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Under the microscope, I found that 
snowflakes were miracles of beauty; and it 
seemed a shame that this beauty should not 
be seen and appreciated by others. Every 
crystal was a masterpiece of design, and no 
one design was ever repeated. When a 
snowflake melted, that design was forever 
lost. Just that much beauty was gone, without 
leaving any record behind. 


Wilson Bentley (1925) 


21.12.21. 
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Preface 


Permanent and seasonal snow covers large areas of terrestrial surface. Unlike liquid 
water, snowpack is a highly reflective surface in the visible region of electromag- 
netic spectrum. The snowpack influences a variety of physical, chemical, biological, 
hydrological, and geological processes. Changes of snow extent and albedo influence 
climate on our planet. Therefore, it is of great importance to study snow properties 
using ground, airborne, and satellite instrumentation. 

This book is an introduction to snow optics. The focus is on the studies of radiative 
transfer in snowpack and the solution of inverse radiative transfer problems related 
to the determination of snow properties from reflected solar radiation. 

My research in snow optics has been inspired by Eleonora Zege, whom I am 
indebted for the introduction to this research filed at the end of last century. My snow 
optics and remote sensing research greatly benefited from cooperation with many 
colleagues around the world including John Burrows, Vladimir Rozanov, Leonid 
Dombrovsky, Wolfgang von Hoyningen-Huene, Teruyuki Nakajima, Teruo Aoki, 
Masahiro Hori, Georg Heygster, Takashi Nakajima, Jason Box, Ghislain Picard, 
Marie Dumont, Maxim Lamare, Florent Domine, Baptiste Vandecrux, Carsten 
Brockmann, Olaf Danne, Simon Gascoin, Marco Tedesco, Alexei Lyapustin, Sergey 
Korkin, Petri Raisanen, Jouni Peltoniemi, Andreas Macke, Stephen Warren, Richard 
Brandt, Tom Grenfell, Don Perovich, Odelle Hadley, Harendra Negi, Iosif Katsev, 
Alexander Prikhach, ^ jio Di Mauro, and Michael Mishchenko. 

While writing this book, my research was supported by the European Space 
Agency and Japan Aerospace Exploration Agency. 

Ithank my family for having endured my preoccupation with working on the book 
over the past several years. 


Darmstadt, Germany Alexander Kokhanovsky 
June 2021 
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Chapter 1 A) 
Microphysics and Geometry of Snowpack ш 


1.1 Ice Grains in Snow: Size, Density and Shape 


Snow is composed of close-packed ice crystals suspended in air. Most of crystals 
in snow have irregular shapes. Although falling snow crystals often have a striking 
symmetry (Bentley and Humphreys 1962). There are multiple ways to classify the 
dimension of crystals. Most common approach is to measure and report the largest 
dimension d of a crystal. Depending on the value of d, a given snowpack can be 
classified in six categories given in Table 1.1 The first three categories are most 
common. Snow crystals originate from atmosphere. Very fine size of crystals may 
serve as an indication of a recent snowfall because crystals grow with time. Snow 
metamorphism is driven by gradients in vapor pressure, which in turn are driven by 
temperature gradients. Small temperature gradients (less than 10 degree per meter) 
result in small vapor pressure gradients and slow grain growth within the snowpack. 
With time, crystals become more rounded. This is because the vapor diffusion within 
the snowpack causes a loss of mass from points on individual snow grains to gains 
in mass in hollows. The large temperature gradient induces a large gradient in vapor 
pressure, such that water vapor moves from warmer areas of the snowpack with 
relatively higher vapor pressures across pore spaces to colder areas of the snow- 
pack with lower vapor pressure. These conditions produce angular or faceted grains, 
which may later develop steps and striations on their surface, resulting in cup-shaped 
crystals with a hollow centre that generally range in size from 3 to 8 mm. Under very 
favourable conditions, individual grains can be larger than 15 mm. The ice crystals 
in snow can be solid, hollow, broken, abraded, partly melted, rounded or angular. 
The surface of crystals can be rimed, stepped or striated. The rounded facets can 
be present as well. The crystals can be bounded, unbounded, clustered, or arranged 
in columns. The ice crystals efsnow can be randomly distributed or partially (say, 
in horizontal plane) oriented. The diversity of crystals in snow is demonstrated in 
Fig. 1.1a, b. 

Each crystal in snowpack differs from other crystals both in terms of shape and 
size. Therefore, the precise characterization of the snow microstructure is much 
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Table 1.1 The international classification of snow grain sizes (the largest dimension of ice grains) 
(Fierz et al. 2009) 























Type Term Predominant maximal dimension of crystals (mm) 
1 Very fine <0.2 

2 Fine 0.2-0.5 

3 Medium 0.5-1.0 

4 Coarse 1.0-2.0 

5 Very coarse 2.0-5.0 

6 Extreme >5.0 








more complex as compared, e.g., to the case of water clouds composed of spherical 
droplets. In case of spherical polydispersions one can introduce the average radius 
of particles 


a= | “ወመ (1.1) 


where a is the radius of particles, f (a) is the particle size distribution normalized as 
follows: 


J f(a)da = 1. (1.2) 
0 


In earlier works related to snow optics, it has been a standard approach to model 
the snow optical properties assuming the spherical shape of particles. This made it 
possible to simplify modelling and derive useful theoretical results for snow radiative 
properties on the basis of well-known results of Maxwell theory as applied to spher- 
ical scatterers. The snow microstructure has been represented by a single parameter 
called the effective grain radius: 


Јаз уба) аа 

as = —— (1.3) 
f а? f (a) da 
0 


or the effective diameter de = 2a,r. The radius aep can be derived from spectral 
reflectance measurements of snow cover. Therefore, it is often referenced as an 
optical grain radius. 

One must be aware of pitfalls of such a simplified approach for snow optics 
problems. Indeed, large transparent ice spheres are characterized by rainbow and 
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Fig. 1.1 а Fresh snow crystals. Courtesy: D. McCoig, b Wind-packed snow from Alaska’s North 
slope (USA). Courtesy: Electron and соо со mper ешш Кезеагсһ 
Service, U. S. Department of Agriculture (ht ነ 5- 
of- a ys als/) 
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glory effects similar to those observed in water clouds. However, such effects have 
not been observed for natural snow covers. Therefore, the account for the irreg- 
ular shape of particles is essential for correct calculations of spectral, angular and 
polarization characteristics of light reflected and transmitted by snow layers. Mete- 
orological classification of snow crystals has been performed by Magano and Lee 
(1966). Ishimoto et al. (2018) have derived sizes and shapes of crystals from X-ray 
computed microtomography imagery. Snow optical properties are modelled using 
various simplified shapes of ice crystals such as hexagonal plates and columns, Koch 
fractals, Voronoi crystals, dendrites, etc. The rough crystals are considered as well. 

Yet another important parameter of a snowpack is the volumetric concentration c 
of ice crystals defined as 


<|= 


(1.4) 


Here N is the number concentration of ice crystals in unit volume of snow 
(measured in cm ^?) and V is the average volume of grains. It follows for spherical 
particles: 


оо 


V= [vereda (1.5) 


0 


where v(a) = 4л a? /3 is the volume of spherical particles with the radius a. In case 
of nonspherical particles, one can use the following equation to derive V: 


оо 


V= “ጋው (1.6) 


0 


where f (v) is the volume distribution function. The snow density p; is related to the 
volume concentration of snow grains c: 


Ps = Сб, (1.7) 


where p; is the density of bulk ice. The value of p; is equal to 0.917 естт at 0 °C and 
standard atmospheric pressure. The snow density is often measured in the filed with 
typical values in the range 0.1--0.4 gem™. The value of volumetric concentration c 
is often in the range 0.1—0.4 and snow porosity р = /—c is in the range 0.6-0.9. 
For such high concentrations of ice grains one can observe some order and corre- 
lation in the position of ice grains. This in principle may influence the snow optical 
properties. However, corresponding effects are usually ignored (at least in the visible 
and near infrared parts of the electromagnetic spectrum). The snow density can be 


1.1 Ice Grains in Snow: Size, Density and Shape 5 


derived weighting a given volume of snowpack. Therefore, the procedure is less 
involved as compared to the complete characterization of snow microstructure. 
1.3 TheSnow Specific Surface Area 


The specific surface area (SSA) o of any solid including ice is defined as total surface 
area S of a material per unit mass M. It follows for snow: 


eS css (1.8) 


Let us assume that snow grains are not in direct contact. Then one derives: 


5 = NS, М = No. V. (1.9) 
Here, 
oo 
s= f stes (1.10) 
0 


is the average surface area of ice grains, f (s) is the surface area distribution function. 
Therefore, it follows: 


o= . (1.11) 
piV 
If grains touch each other and overlap, the accuracy of equation for the SSA 
presented above decreases. The theoretical calculation using Eq. (1.11) provides an 
upper value of SSA as compared to in situ specific snow surface area measurements 
(say, using the methane adsoption technique). The measurements of SSA can be used 
to estimate the effective diameter of grains defined as: 





6V 
ар = —. (1.12) 
5 
Namely, it follows: 
6 
def = ; (1.13) 
o pi 


This formula makes it possible to derive the effective grain size for a given snow- 
pack using snow specific surface area measurements. On the other hand, knowing 
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Fresh new snow 50-70 

Damp new snow 100-200 
Settled snow 200—300 
Depth hoar 100-300 
Wind packed snow 350—400 
Firn (granular) 400-830 
Very-wet 700—800 
Glacier ісе 830-917 








the effective grain size (say, from optical measurements), one may derive SSA: 


Ü 


с = —, 
der 


(1.14) 


where ¢ = 6/p;. One shall recall that formulations presented here have larger uncer- 
tainties for snow having larger values of snow density due to the presence of many 
touching and overlapping grains. The typical values of SSA for different types of 
snow are presented in Table 1.2. 


13 TheSnow Water Equivalent 


The snow water equivalent (SWE) is defined as the depth of water that would theo- 
retically result if the entire snowpack is melted instantaneously. It can be estimated 
as follows: 


SWE = kh. (1.15) 


Here h is the depth of the snowpack, к = p0;/py, Ps is the snow density and 
оь = 1 g/cm? is the water density. The SWE is more relevant for hydrological 
applications as compared to the snow depth because the snow of the same depth can 
contain different amount of water depending on its density. In particular, it contains 
less amount of water at the start of winter season as compared to the snow of the 
same depth at the end of winter season due to t ow densification processes. 

In remote sensing applications the value Кө! usually unknown in advance. 
Therefore, it is parameterized as the function of the snow depth A. In particular, one 
may use the following approximation (Sturm et al. 2010): 


K — Kmin + (Kmax Е Kmin)(l = e ^), (1.16) 
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Table 1 e minimal and г zi 
Е ns Snow class Kmin Kmax y,cm 
maximal densities and 
values of the coefficient y for Alpine 0.2237 0.5975 0.0012 
different snow types (Sturm Maritime 0.2578 0.5979 0.0010 
et al. 2010) = 
Prairie 0.2332 0.5941 0.0010 
Tundra 0.2425 0.3630 0.0049 
Taiga 0.2170 0.2170 0.0000 














Here y is the fitting parameter, Kmin and Kmax are minimal (at the start of the winter 
season) and maximal values of the coefficient x. Typical values of these parameters 
are given in Table 1.3. More accurate formulation of the dependence « (h) includes 
the effects of snow aging, wind, and temperature. In particular, snow aging can 
be accounted for multiplying the exponent in Eq. (1.16) by exp(— Vt), where Ҹ is 
the fitting constant and f is time. 

The snow depth, the primary source of spatial variability of snow water equivalent, 
is an important parameter for many practical applications including various opera- 
tions in snow-covered terrain and also for hydrological models. It is usually measured 
by the snow rules over snow boards, if available. The ground—based tripoid lidar 
systems (terrestrial laser scanner, TLS) are used to monitor the snow depth with a 
high accuracy (5 cm) in automatic way. The snow depth is found via subtraction 
of snow—free from snow-covered surface datasets (Deems et al. 2013). Many TLS 
systems operate at 1550 nm, while most airborne lidar systems work at 503 and 
1064 nm. The systems operating at longer wavelengths have generally smaller pene- 
tration depths, and, therefore, smaller errors in the snow depth estimation related 
to laser light scattering inside snow. On the other hand, small reflectivity of snow 
at near—IR wavelengths limits the operation range. Typical TLS systems deployed 
to characterize snow and ice are significantly range limited («150 m) due to very 
low snow reflectance at 1500 nm. One can use the wavelength 1064 nm to increase 
the snow reflectivity and operation range. However, TLS systems with a 1064 nm 
wavelength are not safe at close range using common pulse and scan rates (Deems 
et al. 2013). 

Both snow depth and SWE are important characteristics of snowpack needed 
for numerous applications. Therefore, they are measured routinely in field and also 
derived using remote sensing techniques. Sample time series of SWE and snow depth 
are given in Fig. 1.2. 


1.4 Layered Nature and Complex Geometry of Snow Fields 


Snow belongs to the class of vertically inhomogeneous turbid media. Its microphys- 
ical (e.g., snow grain size, shape, orientation of grains and density) properties and 
also pollution load change with distance from the snow top. The layered nature of 
snow is due to occurrence of subsequent snowfalls and possible dust falls and other 
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Fig. 1.2 Sample time series (2000-2002) of SWE and snow depth from the Rex River (WA, USA) 
SNOTEL station. Observations of / at times when SWE is zero are spurious (from Hill et al. 2019) 


pollution events. The measured snow grain size, shape, density, and temperature 
profiles for 1.0--1.6 m snowpack are shown in Fig. 1.3. A series of melt/freeze and 
dust loaded layers and also presence of wetting fronts are clearly seen. Melt/freeze of 
ice layers are consistently observed beneath dust layers. One can see that density can 
increase or decrease with snow depth/location in snow layer depending on snowpack. 
The vertical profile of grain size shows large fluctuations with general tendency of 
smaller snow crystals to be located closer to the snow top. 

In addition, the snow surface is characterized by the presence of sastrugi, snow 
dunes, etc. Under the action of steady wind, free snow particles can accumulate and 
drift like the sand grains in barchan dunes. Sastrugi are sharp irregular grooves or 
ridges, which are formed on a snow surface by wind erosion, saltation and deposition 
of snow particles. Sastrugi can have a size of several meters and height 20-30 cm 
(sometimes 1.5 m) with often presence of sharp ice crystals at the top. They are 
oriented parallel to the prevailing wind direction. Sastrugi are similar to water waves 
on the surface of ocean. The example of sastrugi is shown in Fig. 1.4, where snow 
shadowing and brightening effects caused by sastrugi are illustrated. These features 
are clearly seen in satellite imagery. They influence albedo of snow fields and have 
important climatic effects. 
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Fig. 1.3 Grain size, density, and temperature profiles with visual stratigraphy identifying prominent 
features. Pits are labeled with corresponding pit number and date of collection. Pictures from these 
field days are shown below. The buried polluted snow layers are clearly seen (from Skiles and 
Painter 2016) 
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Fig. 1.4 Sastrugi in Antarctica (from NOAA Photo Library—anta0183, Public Domain, https:// 
commons.wikimedia.org/w/index.php?curid-66437699). Diamond dust above sastrugi and also 
sparkling crystals on the surface of sastrugi are present 


1.5 Snow Impurities: Soot, Dust, and Algae 


Natural snow can contain liquid water (wet snow) and various impurities including 
living matter (see Figs. 1.5, 1.6 and 1.7). Snow absorbs certain polycyclic aromatic 
hydrocarbons which are organic pollutants known to be toxic and carcinogenic 
(Nazarenko et al. 2017). The impurity load is characterized by the weight of impu- 
rities divided by the weight of a snow sample. The unit ng/g can be substituted by 
ppbw (parts per billion (107°) weight). The units ug/g or ppmw (parts per million 
(1075) weight) and mg/g or pptw (parts per thousand (10^?) weight) are often used 
for heavily polluted snow. Di Mauro et al. (2015) reported the values of dust concen- 
tration in snow in the range 1—325 ug/g in European Alps after dust deposition 
events. This leads to reddish snow observed over Alps after dust falls originated in 
Sahara. Painter et al. (2012) have reported very high (0.2--5 mg/g) annual end-of- 
melt season dust concentrations at the selected subalpine and alpine sites for time 
period 2005—2010. The sites are shown in Fig. 1.5. The snow pollution is clearly 
seen by a naked eye. It is well known that not only dust but also soot present in snow 
leads to climate forcing via snow albedo (Hansen and Nazarenko 2004; Hadley and 
Kirchstetter 2012). 

The colored (red, green, etc.) algae can present in snow. The examples of snow 
algae are shown in Figs. 1.6 and 1.7. The characteristic algae absorption peaks are 
summarized in Table 1.4. Cook et al. (201778, b) have proposed the physical model, 
which predicts solar light spectral reflectance and albedo of snowpacks contaminated 


1.5 Snow Impurities: Soot, Dust, and Algae 11 


a) 





Fig. 1.5 a Senator Beck Study Plot (alpine site), 12 May 2009 and b Swamp Angel Plot (subalpne 
site), 13 May 2009. Energy balance/radiation towers are depicted. The total dust loading was 
54.6 g/m? for year 2009 (from Painter et al. 2012) 


12 1 Microphysics and Geometry of Snowpack 





Fig. 1.6 Green snow algae (Gray et al. 2020) 


with variable concentrations of red snow algae with varying diameters and pigment 
concentrations. They also estimated the influence of algae on snowmelt. 
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Fig. 1.7 Red snow algae on Harding Icefield in Alaska (Segawa et al. 2018). Also microscopic 
view of the dominated snow algae is given on the right upper panel (scale bar is 20 microns). Segawa 
et al. (2018) have performed single—cell polymerase chain reaction (an isolated cell is shown in 
right lower panel) for the characterization of collected snow algae 




















Table 1.4 Spectral absorption features of algal pigments (Cook et al. 7a) 
Pigment Absorption peaks, nm | Comments 
Chlorophyll a 440,680 Two narrow absorption bands 
Chlorophyll b 475,650 Two narrow absorption bands 
Primary caratenoids | 480 One broad band 
Secondary carotenoids | 460 One broad band 
Phycocyanin 610 One broad band (350-700 nm) 
Phycoerythrin 450, 525, 575 “Table-shaped” absorption spectrum with a 


sharp increase in absorption at 420 nm, plateau 
with a characteristic “triple-peak” morphology 
to a sharp drop at 580 nm (characteristic for red 
algae) 








It is useful to distinguish two classes of snow surfaces: Case-1 snow (fresh not 
polluted snow) and Case-2 snow (polluted snow). This is of importance both for 
the theoretical modelling and for the determination of snow properties from optical 
measurements. The peculiarities of the spectral reflectance and color of polluted snow 
is largely determined by various contaminants (see Figs. 1.5-1.7). The reflection of 
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light from clean snow is determined by the ice grain size/shape distributions and 
optical constants of bulk ice. 


1.6 Optical Constants of Ice 


The spectral reflectance of polluted snow is determined by the size/shape distributions 
of ice grains, local optical properties of impurities and complex refractive index 
m = n — ix of the bulk ice. The spectral behavior of optical constants (n, x) of 
ice is given in Fig. 1.8. One can see that the value of n slightly decreases with the 
wavelength. It changes from 1.325 at the wavelength À = 0.35 um to 1.227 at እ 
= 2.5 jum. The variation of the imaginary part of ice refractive index x covers six 
orders of magnitude. The value of x is very low in the visible. This leads to the fact 
that the absorption of light by clean snow in the visible (as far as light reflectance is 
of concern) can be neglected and snow appears white for human eyes. Such a high 
reflectivity of snow in the visible is of importance for climate change studies (Hansen 
and Nazarenko 2004). The snow reflectance is substantially reduced in the infrared 
region of the electromagnetic spectrum due to larger absorption of light by ice grains 
in this case. 

Water has no absorption bands centred in the visible or near-UV, 0.2-0.7 um. 
The observed absorption by water and ice in this region is the tail of near-IR vibra- 
tional absorptions, decreasing with decreasing wavelength to extremely small values 
(Warren 2019). In weakly absorbing regions, the absorption coefficient о = 4л x /À 
is measured by transmission (Grenfell and Perovich 1981). Such measurements lead 
to large uncertainties at small vales of x at À below 0.6 jum. The results presented in 
Fig. 1.8 contain data presented by Warren and Brandt (2008) with modifications for 
the value of x (A) in the spectral range 0.35--0.6 jum proposed by Picard et al. (2016). 
There is still some uncertainty in the value of x (А) below 0.6 рт. The value of 
X (A) is essentially zero in the visible for some purposes. However, small absorption 
of light by snow in the visible matters for computation of photochemical fluxes in 
snow. Also there is some impact (below 1%) on the snow albedo in the visible. The 
imaginary part of ice refractive index is shown in yet another scale in Fig. 1.9. One 
can see two absorption peaks at 1.5 and 2.0 jum. The horizontal line indicates the 
spectral region (below ~1.15 jm), where the absorption of light by snow is rela- 
tively weak as far as studies of snow local optical properties are of concern. The 
bands located at 1.5 and 2.0 jum include overlapping overtones and combinations 
of the fundamental absorption modes. These vibrational modes, because they are 
intramolecular, are also seen in liquid water and water vapor, shifted somewhat in 
frequency. The corresponding shifts (also at shorter wavelengths, see Fig. 1.9) can 
be used to identify the abundance of various thermodynamics states of water in snow 
(Green et al. 2006). 

Ordinary ice is birefringent. Therefore, light refraction depends on the angle 
between the c-axis of the crystal and the direction of propagation and polarization of 
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Fig. 1.8 Spectral complex refractive index of ordinary hexagonal ice Ih (Warren and Brand 2008; 
Picard et al. 2016) at temperature 266 K 
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Fig. 1.9 The imaginary part of ice refractive index (Warren and Brand 2008) 


incident light. However, the difference in refractive indices of ordinary and extraor- 
dinary rays is very small (0.196) and ice birefringence can be often ignored in snow 
optics studies. One can neglect the temperature dependence of the complex refractive 
index of ice in the UV, visible and near IR parts of electromagnetic spectrum. This 
is not the case in the mid-IR and far-IR bands. 
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Chapter 2 A) 
Local Optical Properties of Snowpack get 


2.1 Geometrical Optics of Large Spherical Particles 


2.1.1 Light Scattering 


Snow crystals are much larger as compared to the wavelength А. of incident light. 
Therefore, the interaction of light with a single crystal can be studied in the framework 
of geometrical optics. To simplify, we consider at first main results of geometrical 
optics as applied to spherical scatterers with the radius a and the relative complex 
refractive index m — n—ix. Let us assume that x — oo and 2x(n — 1) — oo, 
where x = ka is the size parameter and k = 2л. This makes it possible to ignore 
phase shifts and possible interference effects. 

The interaction of light beam with a sphere with the radius a much larger as 
compared to the wavelength A of incident light is demonstrated in Fig. 2.1. The 
incidence angle i shown in Fig. 2.1 varies from 0 (incidence along the diameter of 
sphere) to 90? for the grazing incidence. For a given incidence angle, the reflected 
and refracted beams originate. The reflected beam leaves the particle at the reflection 
angle equal to the incident angle i. The refracted beam propagates in the direction 
specified by the refraction angle r as it follows from the Snell's law (sin r — 
n sin i, where it is assumed that ዘ (see Appendix)). One can see that the path 
of the ray becomes longer as compared to the undisturbed ray propagating in the 
direction specified by the angle i. The ray deviates in the direction of the center of 
a particle (focusing effect of the particles). This means that possible impurities (say, 
soot) located inside of scatterers act as more effective light absorbers (due to larger 
absorption on longer paths inside particles) as compared to the same impurities 
outside scatterers. Therefore, external and internal mixture of impurities in snow 
will lead to different light absorption in snow and also differences in snowmelt 
processes and timing. The refracted light is partially refracted for the second time 
and leaves the particle at the angle i. Another portion of the ray is reflected inside 
the particle and the process is repeating infinite number of times forming a complex 
scattering pattern around the particle. In case of absorbing particles, one needs to 
© Springer Nature Switzerland AG 2021 19 
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Fig. 2. Ray tracing ina 
sphere 





И, 


account for the absorption processes on the path / of the ray inside the particle. 
This is done taking into account the exponential attenuation of the ray intensity J on 
the distance L according to the following law: J exp(—o L), where о = NE is the 
bulk ice absorption coefficient, x is the imaginary part of ice refractive index at the 
wavelength A. 

Let us consider first the reflected light beam and introduce the complementary 
angles 


us х= x тыр, Я 


The use of these angles makes the final equations somewhat simpler. In particular, 
it follows from the Snell's law at x «& n: 








cost! = ČZ, (2.2) 
n 
This also means that 
m cost? 
sini 1-[ |. Q.3) 
n 


A finite pencil of light is characterized by the spread in the azimuthal angle dọ and 
also dt. The flux of energy in this pencil is /oa?costsinvdtdq for a light coming in 
the direction specified by the angle x (Shifrin, 1951; van de Hulst, 1957). Here Io is the 
intensity of incident light. The reflected part of the energy is Jo R ja? costsinzd tdo; 
for the light polarized perpendicular (j = 1) or perpendicular (j = 2) to the incidence 
plane. It follows (Fresnel 1831; van de Hulst 1957): 


š . 12 
SINT — MSINT 


(2.4) 


msint + sint’ 


. А 2 
msint — sint! 
sint + msint’ 











It is known that x «& n for ice in the visible and near infrared. Therefore, one can 
assume that m = n in Eq. (2.4). 

The emergent pencil of light spreads into a solid angle dQ = sin@d@dq. The 
intensity of reflected light beam at a large distance r in the direction specified by 
the scattering angle 0 is found dividing reflected light flux by the area dA = r?d ፔ 
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Namely it follows: 
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ዞራ sin2t Q.6) 
^ 2usin0 ` 
and 
ወፀ d0' 
v= — =|—}. (2.7) 
dt dt 
Here, the angle 
0! = 2т (2.8) 


describes the deviation of the reflected ray from the initial direction. It defines the 
scattering angle = аве (0-180ግ: 


ፀ= 2лК + €] (2.9) 


where K is an integer and = —1 or q = 1. It is easy to show that for other beams 
(two times refracted, refracted-reflected-refracted, etc.), it follows: 


1,.; = Io(1— Ri) R Dar ሞኙ L9 as. (2.10) 
[3. case of absorbing spheres and x <ፎ n, it follows: 


1; = b(1 — Rj) RT! D(a/r)2 exp(—al), (2.11) 





where / = zsin(r’), а = 2kx, k = 2л, 2 = 2as,s=1,2,3,.... 
Summing up, the total geometrical optics scattered miens) be written as 


оо 
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where 
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ipea e D, (2.13) 
= (1 — Rj) Rt exp(—cssint’, (2.14) 
x = Ка, c = 2aa. (2.15) 


The geometrical optics scattering cross section is defined as an integral of intensity 
given by Eq. (2. "CET it follows: 


2л л 


1 a. 
Са, = ú de | I*(0) r^sin0d0. (2.16) 
0 
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0 


It follows from Eqs. (2.16) and (2.12): 
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j=l 
where 
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р, = У (1 — Rj) RT exp(—cs sint’). (2.18) 
5=1 
We derive performing ће summation in Eq. (2.18): 
ን —csint' 
(1 Е к.) е 
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Therefore, it follows (Kokhanovsky and Zege 1995): 
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2 1- R.) e cint 
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Опе can see that the geometrical optics scattering cross section of an absorbing 
spherical particle can be expressed via a single integral over the incidence angle т. 

As a matter of fact, scattered light around е spherical particle can be presented 
as a sum of contribution due to two di processes: diffraction and reflec- 
tion/refraction of rays at the surface of particles. The cross section shown in Eq. ( 
accounts just for the reflection and refraction processes. The intensity of diffr: 
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light can be presented by the following equation (van de Hulst, 1957): 


Ioa? J? (0x) 
а _ 1 
Ica =, 7 022° (2.21) 
where J, (0x) 15 the Bessel function. It follows from Eqs. (2.16) and (2.21): 
oo 
C. = 2ла? | Ji Q)y "ወ (2.22) 
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where у = Өх and we have assumed that x — oo and sin0 ~ 0 at small scattering 
angles. We extended the upper limit integration to infinity taking into account that 
the scattering angle is fixed and small but x — oo. The integral in Eq. (2.22) can be 
evaluated analytically: 


oo 
2 | доу = 1. (2.23) 
0 
Then it follows: 
C“ =a’. (2.24) 


Therefore, total scattering cross section is given by the following expression: 


Csca = та?(1 + W), (2.25) 


where 
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Let us introduce the scattering efficiency factor: 
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It follows for this factor: 
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1 
Оса = 1+ 2 | (шу + w2)sin2tdt. (2.29) 
0 


The scattering efficiency factor can be also derived using Maxwell electromag- 
netic theory (Maxwell 1873; Mie 1908). The derivations are much more involved as 
compared to our discussion valid for the case of large spherical scatterers. The inter- 
comparison of geometrical optics and Mie theory calculations of О, „„ for monodis- 
persed ice spheres with radii 100, 500, and 1000 microns in the spectral range 0.3-2.5 
microns is given in Fig. 2.2. As one might expect, it follows that the accuracy of 
geometrical optics approximation increases with the size parameter of particles. The 
largest deviations occur for 100 um spheres (fine grained snow). The high frequency 
oscillations in Fig. 2.2 are due to interference of various geometrical optics rays (van 
de Hulst, 1957) ignored in the theory described above. They do not appear if one 
takes into account the polydispersity of snow grains. 

One can derive from Eq. (2.27) at c = 0 (nonabsorbing particles): w; = 1 and, 
therefore, W — 1 (see Eq. 2.26). This means that 


Osca =2 (2.30) 


1000um 


scattering efficiency factor 





04 06 08 10 12 14 16 18 20 22 24 
wavelength, um 


Fig. 2.2 The intercomparison of geometrical optics (symbols) and Mie theory (lines) calculations 
of Qsca for monodispersed ice spheres with radii 100, 500, and 1000 microns 
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in this case. This result is very important. It shows that the scattering efficiency factor 
does not depend on the refractive index of particles and their size in the geometrical 
optics limit for large nonabsorbing spheres. Equation (2.30) is valid for nonspherical 
particles as well (van de Hulst 1957). 

Let us use the following approximation valid at small values of the attenuation 
parameter c: 


wj = 1—csint’. (2.31) 
Therefore, it follows: 
Qsca = 2 — Ac, (2.32) 
where the integral 
z 
A=  “ዘሠሥ (2.33) 
0 


can be evaluated analytically using Eq. (2.3) and the substitution: и = cos*t, du = 
sin2tdt. The answer is: 


2 EX 
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One can see that A — 2 аѕ п — 1. Equation (2.29) can be parameterized as 
follows: 


Qsca = 1+ p + (1 — p)exp(—ve), (2.35) 


where v is the fitting parameter and the integral 
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p= 2 (Ri + R;)sin2rdc (2.36) 
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can be evaluated analytically at x ««n. Namely, it follows: 
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where p; = (—1, —1, —3, 7, —9, —13, —7, 3). 

The fitting parameter v can be found from the condition that Eq. (2.35) must 
coincide with Eq. (2.32) as c — 0. It provides the correct behavior of Csca at small 
values of c. Then it follows: 





መው AN (2.38) 


where A is given by Eq. (2.34). The value of р given by Eq. (2.37) can be 
parameterized as follows at n = 1.2--1.4: 


p = 0.0123 + 0.1622(n — 1). (2.39) 


The intercomparison of calculations according to Eqs. (2.29) and (2.35) is shown 
in Fig. 2.3. It follows that analytical Eq. (2.35) can be used with the error less than 1% 
at n = 1.31. One concludes from Eq. (2.29) and Fig. 2.3 that the scattering efficiency 
factor for large spherical scatters changes from 2 for nonabsorbing particles to 1 + 
p for strongly absorbing particles. The value of o depends on the refractive index n 
and close to 0.06 at n = 1.31. The same is true for large nonspherical particles. 

The normalized scattering light intensity 


— — Eqs. (2.35), (2.37) 
geometrical optics 


scattering efficiency factor 





attenuation parameter 


Fig. 2.3 The dependence of the scattering efficiency factor on the attenuation parameter c calculated 
using Eqs. (2.35), (2.37) (red line) and geometrical optics Eq. (2.29) (black line) at n = 1.31 
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+t 
p*(9) TI. 2: (2.40) 


is called phase function. It follows: 
оо 
1 = » is (2.41) 
s=0 


and ij, for each о ел by Eq. 81... of numerical calculations of ће 
phase function for large spherical particles are given by Shifrin (1951) and van de 
Hulst (1957). Zhou et al (2003) provide the code, which can be used for the calcula- 
tions of phase functions of large spherical particles in the framework of geometrical 
optics. 

The total phase function, which accounts for both geometrical optics and 
diffraction components, can be presented in the following way: 


Са p^ (0) + Chap! (0) 








0)= ; 2.42 
р(0) сї Се, (2.42) 
where 
d x?F(0) 
p° (0) = J (2.43) 
and 
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F(0) = E m (2.44) 
The phase function p(0) is normalized as follows: 
Lf M 
; | Р®)зїпөдө = 1. (2.45) 
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Similar equations hold for р“: (0). The phase function gives the conditional prob- 
ability of light scattering in a given direction and can be used to calculate various 
statistical parameters relevant to the p ation of light in turbid media. In particular, 
the average cosine of scattering angl n be derived: 
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g= ; | p@sinocossae. (2.46) 
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Let us derive useful approximation for the average cosine of scattering angle in the 
framework of geometrical optics. Then it follows (van de Hulst 1957; Kokhanovsky 
and Zege 1995): 








1+ We? 
Sede (2.47) 
l+ W 
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The inter-comparison of asymmetry parameters calculated using Eqs. (2.47-2.49) 
and Mie theory is shown in Fig. 2.4. One can see that size dependence of the 
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Fig. 2.4 The intercomparison of asymmetry parameters calculated using Eqs. (2.47-2.49) 
(symbols) and Mie theory (lines) for the radii of particles equal to 100, 500, and 1000 m 
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asymmetry parameter is weak in the visible. The results of geometrical optics and 
Mie calculations almost coincide for large values of size parameter. Therefore, fast 
geometrical optics calculations can be used instead of tedious Mie computations for 
large scatterers. 

It follows at c = 0: 


(1— Rj) cos2(t = ፻፲) + 2፳7፡ዐ፡27(1 = cos2t') 
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In the opposite case of strongly absorbing particles (ር — oo) one derives: 
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(2.56) 


where q; — (—3, 13, —89, 151, 186, 138, —282, 22, 25, 25, 3, 3). 
Unfortunately, analytical integration for the asymmetry parameter g at any c is 
not possible. However, one can use the following uniform approximation: 


g = 8 — (፪% — 80) exp(—cn), (2.57) 


where go; = g(oo), go = g(0) and the value of 7 can be estimated, e.g., using the 
asymptotic behavior of g given by Eq. 4 at small values of c. Then it follows: 


g = go + cX, (2.58) 
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X = P + Т, (2.59) 
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М; = 1 — 2Rjcos2t' + Rš. It follows from Eq. (2.57) at small values of c: 





g = go + ርዛ(፪% — 80). (2.62) 
Therefore, one derives from Eqs. (2.58) and (2.62): 
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where X is defined by Eq. (2.59). The parameters go, ድሬ and n depend only on the 
real part of ice refractive index at x «& n. They can be parameterized as follows: 


go = 1.006 — 0.3641 (n — D), вә = 1.008 — 0.1(n — 1), 
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Fig. 2.5 The dependence of the asymmetry parameter of ice spheres on the attenuation parameter 
at the real part of ice refractive index 1.31 and radius of particles equal to 1000 um. The wavelength 
is equal to 1 um. The results of calculations using the geometrical optics approach (Eq. 2.47), 
uniform approximation (Eq. 2.57), and linear approximation (Eq. 2.62) are shown 


п = 0.3639 + 1.676(n — 1) — 1.6284(n — 1)? (2.64) 


at n = 1.2-1.4. The accuracy of Eq. (2.57) is presented in Fig. 2.5. The error is below 
0.596. The accuracy of linear approximation (see Eq. 2.62) is also shown. It can be 
used for the attenuation parameters c smaller than 0.5. 


2.1.2 Light Absorption and Extinction 


The intensity of transmitted light for the j-th polarization interacted—s-time—with 
the-surface-of-the-particle is given by Eq. (2.11). In particular, it follows for the 
transmitted beam at s = /: 


I; = 601 — Rj)’ D(a/rY exp(—csint’). (2.65) 
It follows in the absence of absorption: 


Ij = ከ0 - RjYD(ajrY. (2.66) 
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Therefore, the absorbed light intensity on the ray path can be presented as 


1 (s = ]) = 1(1 ШЕ Rj)’ D(a/ry (1 — exp(—csint’)) (2.67) 


at s = 1. It follows at s = 2: 
14% (s = 2) = Ip(1 Е) (а/г), (1 — exp(—2csint’)). (2.68) 
One derives summing all rays: 


In — R )Da2(1 — exp(—csint’)) 











pv 2.69 
J (1 — Rjexp(—csint’))r? VPE 
The absorption cross section is defined as 
2л л 
1 
Cu c а | аф J Q^ (0) + 12^ (0))r2sin0d0. (2.70) 
02 0 
Therefore, one derives: 
Caps = Mra’, (2.71) 
where 
M fe + тә)зїп2та Euer Q.72) 
= — | (m, +m))sin2tdt, m; = TT . . 
2 ! 2 1 ሸመ Rje cnt 
0 
The value of M coincides with the absorption efficiency factor 
Cabs 
Qars = —5. (2.73) 
ла 
It is easy to show that 
M+W=1. (2.74) 


This means that both absorption and scattering cross section are determined by a 
single geometrical optics integral (e.g., W, see Eq. 2.26). Therefore, we can write: 


Cabs = za? (1 — W), Coca = ла?(1 + W). (2.75) 
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Fig. 2.6 The dependence of the absorption efficiency factor of ice spheres of various radii on the 
wavelength in the near infrared (solid lines-Mie theory, symbols-geometrical optics) 


It follows for the efficiency factors: 
Qavs =1- W, Qsca =1+W. (2.76) 


An important point is that efficiency factors depend on just two parameters in the 
geometrical optics approximation (at x «« n). They are real part of refractive index n 
and attenuation parameter c — od, which is equal to the product of bulk ice absorption 
coefficient and the diameter of particles a. The accuracy of Eq. (2.76) for the 
absorption efficiency factor is shown in i 

The extinction cross section Cext is defined as the sum of scattering Csca and 
absorption Caps cross sections: 


Cext = Csca + Cabs. (2.77) 

Therefore, we arrive to the following equation for the extinction cross section : 
Cag = 2ma2. (2.78) 
This equation is valid at any c in the framework of the geometrical optics approx- 


imation (x — oo, 2x(n — 1) — oo). This result is very important. It shows that the 
extinction cross section of large spherical particles does not depend neither on the 
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extinction efficiency factor 
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Fig. 2.7 The spectral dependence of the extinction efficiency factor of ice spheres derived using 
Mie theory at radii equal to 25 (blue line), 50 (green line), and 100 (red line) microns 


wavelength nor on the refractive index of particles. The extinction efficiency factor 
Qe = Cox is equal to 2 for large spherical particles independently of their actual 
size and chemical composition as far as x — oo and 2x (n — 1) — oo, where x is the 
size parameter. The spectral dependence of the extinction efficiency factor derived 
using Mie theory for monodispersed ice spheres with radii 25, 50, and 100 microns is 
shown in Fig. 2.7. The oscillating feature is due to the interference of waves passing 
the spherical scatterer of a given radius. The spectral oscillations disappear for real- 
istic polydispersed media such as snow. It follows that the error of the asymptotic 
approximation (Qext = 2) is better than 1% at the spectral range 0.3-2.5 jum and 
radii larger than 100 microns typical for snow covers (see Table 1.1). 

Let us assume that particles are weakly absorbing and, therefore, c — 0. Then it 
follows from Eqs. (2.72) and (2.32): 





Cabs = bo V, (2.79) 


where B = 3A/2 or 


В = n2(1 — (1 — n 232). (2.80) 
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We can see that В > 1 85 n > 1 (optically soft particles). It follows from 
Eq. (2.79) that absorption cross section of large spherical particles is proportional to 
the volume of particles = id the bulk absorption coefficient o. This remains true 
for nonspherical particles 45 well. However, the coefficient B differs for particles of 
different shapes. It follows for arbitrarily shaped optically soft large particles: 


Cabs = ወሃ. (2.81) 


The useful approximation for the value of Caps valid at any values of n and c can 
be derived from Eqs. (2.35) and (2.77): 


Cans = za? (1 — р)(1 — exp(—vo)). (2.82) 





The parameter v can-be-foundfromthe-condition-that ር must-coincide-with 





y= 2B (2-83) 
ር 30 = o0) ` 

Qus = Cans (2-84) 
^ sta i 






-It follows 


Qavs = (1 — p)(1 — exp(—ve)); m -6 











p = 0:0123 + 041622 = H 286) 





Qs = $+ р F A= р) expre} 87) 
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Fig. 2.8 a The absorption efficiency factor of ice spheres with radii 100, 500, and 1000 um; 
b The scattering efficiency factor of ice spheres with radii 100, 500, and 1000 um (solid line— 
geometrical optics, dash line—linear approximation, symbols—uniform approximation). The linear 
approximation can be used at the wavelengths smaller than 1.2 m 
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Fig. 2.9 The dependence of the absorption efficiency factor of ice spheres on the attenuation 
parameter c at the real part of ice refractive index 1.31 and radius of particles equal to 1000 
calculated using geometrical optics approximation and also linear and uniform approximations. 
The wavelength is equal to 1 ит. Similar accuracy of approximations holds for the scattering 
efficiency factor (Qsca = 2 — Qabs) 


The accuracy of approximations for the values of Оль; and О, са is given in 
Fig. 2.8. It follows that the error of the uniform approximation is better than 2%. One 
concludes from Ege.85) and Fig. 2.9 that the absorption efficiency factor for large 
spherical scatters а ges from 0 for nonabsorbing particles to 1—p for strongly 
absorbing particles. The same is true for large nonspherical particles. The value of o 
is given by Eq. (2.37) at ሂ<<ክ, which is valid assumption for ice in the visible and 
near-infrared regions of electromagnetic spectrum. The absorption efficiency factor 
Qaps is directly proportional to the attenuation parameter c for weakly absorbing large 
spherical particles: Qaps Ac, where the coefficient A depends on the refractive 
index of particles (see Eq Ei 4). 

The parameters introduced above can be used to calculate other light scattering 
characteristics useful in radiative transfer studies. In particular, the results for the 
single scattering albedo wo [bs and the probability of photon absorption (PPA) 
B=1 ГТ of monodisper: pherical particles are shown in Fig. 2.10a, b. 
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Fig. 2.10 The spectral dependence of single scattering albedc ሯ d probability of photon absorp- 


tion b for monodispersed ice spheres with radii 100, 500, а )00 um (symbols- geometrical 
optics, lines-Mie theory) 
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2.20 Local Optical Properties of Snowpack 


2.2.1 Integral Light Scattering and Absorption 
Characteristics of Large Nonspherical Scatterers 


Particles in snow are much larger as compared to the wavelength of the incident light. 
Also they have nonspherical shapes and respective phase shifts are large. Therefore, 
the extinction cross section can be presented in the following way (van de Hulst 


1957): 
Cext = 2X, -e 


where X is the cross section of the particle in the direction perpendicular to the 
incident light. It follows that extinction cross section of a large nonspherical particle 


coincides with that of a sphere with the radius a = 


Let us assume that nonspherical particles are сопу randomly oriented. Then 
it follows for the average value of X: 


, (2.89) 


tə|@l 


where Ф is the average surface area of the particles. Therefore, it follows: Cex: = 
The phase function, asymmetry parameter, scattering cross section and absorption 
cross section of nonspherical particles larger as compared to the wavelength of the 
incident light can be calculated using the ray tracing approach. This has been done 
for various shapes of particles and corresponding databases are available online. 

The simple geometrical optics approximation for the average absorption cross 
section of randomly oriented nonspherical particles of the same size can be presented 
in the form similar to that for spherical particles: 


Caps = (1 — p)(1 — exp(—we)){), (2.90) 


where c = aD, D = a (27) is given by Eq. (2.89), V is the volume of particles, 


and w is the fitting parameter. One can see that Caps for large strongly absorbing 
nonspherical partilces in random orientation can be presented E] 


Cap; = (d — p){). (2.91) 


It follows in the case of weakly absorbing particles: 
Cabs = Ba V, (2.92) 
where 


B = 5(1—p)w. (2.93) 
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Clearly, one derives: B — 1, о — 0 as n — 1 for particles of arbitrary shapes. 
This means that w — 2 as п > |. It is of importance that parameter B does not 
depend on the size of particles for large weakly absorbing scatterers. The dependence 
of B on the shape of particles has been studied by Libois et al. (2013). The computa- 
tions have been shown that the value of B varies in the range 1.25 (spheres) to 2.09 
(spheroids with aspect ratio 0.5). The following variety of shapes have been consid- 
ered: cylinders, spheres, spheroids, hexagonal plates, fractal particles, and cuboids. 
The analysis of measurements performed by Libois et al. (2013) have shown that 
the most frequent value of B is in the range 1.6--1.9 for natural snowpacks. Taking 
В = 7/4 as a representative value, we derive from Eq. (2.93) at n = 1.31: w = 1.24. 
This constant can be used in conjunction with Eq. (2.90) to study the dependence 
of the absorption cross section of snow grains on the wavelength (under assumption 
that one can neglect the spectral variation of the real part of the ice refractive in the 
spectral range under study). 

The asymmetry parameter for large nonspherical partilces of the same size and 
shape in random orientation can be presented in the way similar to the case of 
spherical scatterers: 





8 = ፪% — (Boo — 80) eXP(—Y c), (2.94) 


where у 15 the fitting parameter, ር = o D. It follows that g = ደሬ at c — oo. The 
value of goo is given by Eq. (2.54) in case of randomly oriented nonspherical convex 
particles. The values of go and y depend on the shape of particles. 

Equations given above are valid for the convex randomly oriented particles of a 
given shape and size. In reality particles of different shapes and sizes present in snow. 
Therefore, Eqs. (2.90), (2.94) must be averaged with respect to the size and shape of 
particles. In particular, it follows for the averaged quantities: 


(Co) = 202), (2.95) 
< Cass >= (1 — p)(1 — exp(-opD))(X), (2.96) 
< g >= go — (Boo — (80)) exp(-oqD), (2.97) 


where (27) is given by Eq. (2.89) for convex particles in random orientation and 


— 6V 
= (2.98) 


is the effective grain diameter. Clearly, one derives: (С) = (Сех) — (Cabs). One 
can assume that (go)is equal to 0.76 as measured in situ for ice clouds in the visible 
(Garret et al. 2001). Generally, the value of (go) increases with the roundness of 
grains. The old snow has larger values of (go) (e.g., around 0.8). The pair (p, q) 
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can be derived from the fitting results derived using geometrical optics approach for 
particles of various shapes and their mixtures. 
It follows for the probability of photon absorption of large nonspherical partilces: 





1 c 
B= 50 — p —expC-apDy. (2.99) 
Snow grains are weakly absorbing in the visible and near infrared. Then it follows: 


< Cas >= Bav, (2.100) 


where the parameter B depends on the shape of particles. 


2.2.2 Integral Light Scattering and Absorption 
Characteristics of Snowpack 


The snow is composed of particles having different shapes, sizes, and orientations. 
Therefore, the snow scattering, absorption, extinction coefficients and asymmetry 
parameter are defined as 





C T 22 ር 
ይ ረረ = N Coens Kats. —NC ku =N Coes ያ= 8 S sca 





| (2.101) 





sca 


Here, N is the number of particles in unit volume of snow and line above the 
symbol means the averaging with respect to the geometrical parameters and orien- 
tation of particles. Equations (2.101) can be used for calculations of other integral 
light scattering and absorption characteristics useful for studies of radiative transfer 
in snow layers. They are listed in Table 2.1. The snow optical thickness (SOT) т is 





Table 2.1 Local optical 
characteristics. The light 
diffusion parameters listed in 
the last three lines are derived Probability of photon | 6 1--ዐወ0 
for weakly absorbing strongly absorption 





Single scattering wo 2с 
albedo 








light scattering media Symmetry parameter | € 1—g 
>1 
(oo ) Transport extinction | ደፁ (1—8) Куса + kabs 
coefficient 





Diffusion exponent k(w > 1) |43(1 — ወ0)(1 — cog) 


Diffusion extinction /'(ወ0 > 1) | Khext 
coefficient 








Diffusion similarity #(ሠ0 > 1) |K/3(1— 2) 
parameter 
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defined as: т = kexth, where h is the snow geometrical thickness and it is assumed 
that kext does not change along the vertical coordinate z. Otherwise, it follows: 


<2 


т = ка, (2.102) 


51 


where 22 gives the position of the upper snow boundary, z; is the vertical coordinate 
of the lower snow boundary and, therefore, h = 22 — zi. 
Introducing the ice grain volume concentration c — N V, we derive: 


kext = = (2.103) 
D 





3 __ 
kabs = UI p)(1 — exp(—apD)) (2.104) 


and Ksca = kext = Карз. 

The volumetric concentration of ice particles in snow is often close to 1/3. Then 
it follows that kext = D Taking into account that the effective diameter of grains 
is usually in the range 0.1-1 mm and kext is in the range 1—10 mm-|, we obtain 
that snow optical thickness is usually larger than 5, where we assumed that the snow 
layer thickness Л > 5 mm. The snow with the depth of 10 cm (or even smaller than 
that for the fine-grained snow) has optical thickness above 100 and, therefore, in a 
good approximation such a snow layer can be considered as a semi-infinite medium 
in the visible range. This thickness is even smaller in the near infrared range, where 
ice also absorbs light. Therefore, photons have smaller probability to survive during 
propagation in snowpack. Such a behavior of snowpack has important implications as 
far as radiative transfer in snow is concerned. Firstly, the direct light beam attenuates 
very fast and can be neglected at h 5 mm or so. Therefore, the internal light field 
in snow originates mostly from አጠ liffuse light, which attenuation coefficient y is 
much smaller as compared to the extinction coefficient kext. Secondly, for most of 
applications snow can be considered as a semi-infinite medium in the optical range. 
Therefore, the influence of underlying surface on snow reflectance can be neglected 
in most of cases. Thin layers of snow (5 mm or so) can be considered as turbid 
layers with large optical thickness. Therefore, simplified asymptotic solutions of the 
radiative transfer equation valid at large value of x (see below) can be used. 

It follows in the visible and near infrared (till roughly 1.24 um for natural snow) as 
discussed above: 


kabs = Bac. (2.105) 
In this range snow is weakly absorbing and main equations simplify. In particular, 


one can neglect the dependence of the asymmetry parameter (g) on the snow grain 
size and refractive index due to small variation of n and also because a D — 0. Then 
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it follows: 
(g) © (до). (2.106) 


Also one derives as а D — 0: 


б 
В = seBD. (2.107) 


к = Ja B(1— go) D, (2.108) 


1 aBD (2.109) 
5 == = = s . 
3Y (1 — go) 


where the parameters В, к, s are defined in Table 2.1. One can also introduce the 
asymptotic flux attenuation coefficient (AFEC) 


y = kkext (2.110) 


and transport extinction coefficient for nonabsorbing snow ( at kabs < < kext): 


ky = (1 — &0)Кехт. (2.111) 


The transport extinction coefficient k, determines the light transmittance through 
finite optically thick nonabsorbing snow layers. It coincides with the radiation pres- 
sure coefficient for nonabsorbing media. This is due to the fact that оу, са is propor- 
tional to the forward carried momentum in the scattered light. It follows for these 
parameters in the weak absorption approximation: 


— 3(1 — 
у= асар, k lue (2.112) 


D 


The parameters mentioned above describe integral light scattering, extinction and 
absorption characteristics of snow layers. The e-folding depth z, = 1/y can be 
presented as 


1 


3c/ a B(1— aD. 


The dependence of ze on the wavelength and size of particles at c = 1/3, go = 3/4, 
B — 1.6 calculated using imaginary part of ice refractive index derived by Picard 
et al. (2016) is shown in Fig. 2.11. The value of z, gives the distance at each the 





Ze = (2.113) 
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Fig.2.11 a The dependence of the e-folding depth on the wavelength and effective grain size for 
clean snow in the visible and near UV. b The dependence of the e-folding depth on the wavelength 
at several values of effective grain size for clean snow in the visible and near UV. c The dependence 
of the e-folding depth on the wavelength at several values of effective grain size for clean snow in 
the near infrared range of the electromagnetic spectrum 
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Fig. 2.11 (continued) 


diffuse light intensity decreases in e-times (in the deep layer regime). This e-folding 


depth can be presented as: z,(A) = 0.454/АР/ усе (À) for the parameters specified 
above. Here, Xice(à) is the imaginary part of ice refractive index (Picard et al. 2016). 


2.2.3 Phase Function 


The snow phase function p(0) describes angular pattern of light scattered by a unit 
volume of snow. It can be calculated assuming a particular shape of particles or 
particular mixtures of different shapes (Liou and Yang 2016). The phase functions 
of spherical particles are characterized by peaks in the rainbow and glory regions, 
which absent for snow. Therefore, the assumption of spherical scatterers can not be 
used for the calculation of snow phase function. 

Various parametrizations of snow phase functions have been used in studies of 
light propagation in snow. In particular, one can use the following parameterization: 


N 
p(0) = У anp, 0) (2.114) 


ዘ=1 
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where 
p, (8) = exp(—b,0) (2.115) 


or 





1-82 
(1 + b? — 2b,cos0)? 


p. (0) = (2.116) 


w 


ә! 


The function (2.116) is called the m © Ween (H-G) phase E P 
other parametrizations or combinations of functions given above can be used. In 
particular, it follows at N = 1 from Eq. (2.114), if one uses the function (2.116): 











1-2? 
p(0) = j5 (2.117) 
(1 + bj — 2bicos0)? 
where b, = g is the asymmetry parameter. It follows at N = 2: 
1-5 1 — bš 
р(60) = Y =+(01- T) =. (2.118) 
(1 + bi ፦ 2bicos0)° (1 + bi ፦ 2b;cos0)? 
The asymmetry parameter for this phase function is 
а= ТЬ + (1 — Y) (2.119) 
and the backscattering fraction 
B= | p(0)sin0d0 (2.120) 
2 
is given by the following formula: 
B=YB,+(-Y)b, (2.121) 
where 
1—bj;| 1+b; 
B = J ችን egi. (2.122) 





2bj 1+5 


The coefficients YT, bı, b; in Eq. (2.118) can be selected in such a way that the 
asymmetry parameter g, the backscattering fraction B and phase function in the 
backward direction 


2.2 Local Optical Properties of Snowpack 47 


a 100 





+ Voronoi phase function 
double H-G phase function 


phase function 


0 20 | 40 60 80 100 120 140 160 180 
scattering angle, degrees 


b 1000 


Fractal particle 
Voronoi particle 


100 


10 


phase function 


0.1 





0 20 40 60 80 100 120 140 160 180 
scattering angle, degrees 


Fig.2.12 a The Voronoi aggregate phase function (size parameter 311, n — 1.31) and double 
Henye-Greenstein phase function at Y = 0.8615, bí = 0.85, b2 = 0.2. b Phase function of ice 
crystals at the wavelength 0.56 um (n = 1.31) calculated using fractal (tetrahedron of second 
generation with the side of 100 wm) (Macke et al. 1996) and Voronoi (the size parameter 311) 
aggregate (Ishimoto et al. 2010) irregularly shaped particles 
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1—Ь, 1— b 
s FUST) j 
Q b) (1 + 20) 





p(w) = Y (2.123) 


coincide with the results derived from the numeric calculations in the framework of 
ray optics approach for particles of various shapes. In particular, the results derived 
for the fractal particles (say, Koch fractals of second generation (Macke et al. 1996)), 
Gaussian random particles (Muinonen 1996), Voronoi (1908) aggregates (Ishimoto 
et al. 2010, 2018; Letu et al. 2016) and other shapes (see, e.g., Libois et al. 2013, 
2014) can be used. 

The intercomparison of the double Henyey-Greenstein (1941) phase function 
and the phase function derived for the Voronoi aggregate is given in Fig. 2.12a. The 
Voronoi aggregate particles are constructed using the spatial Poisson—Voronoi tessel- 
lation to shape modelling of snow grains (Ishimoto et al. 2010, 2018). The discrep- 
ancy of Voronoi aggregate and double H-G phase functions is smaller as compared to 
differences in phase funs for different shapes of particles. In particular, we show 
the intercomparison d. and fractal particles in Fig. 2.12b. They differ consid- 
erably in the angular range 0—30 degrees. However, the behavior is quite similar in 
the backward hemisphere relevant to snow remote sensing. One can use the following 
parameterization of the fractal ice grains phase function, which follows from Eqs. 
(2.114 and 2.115) at bj = 0.087 deg !, b; = 0.014 deg l, a; = 11.1, аз = 1.1: 


p(0) = 11.1exp(—0.0870) + 1.1 exp(—0.0140). (2.124) 


This equation can be used only outside diffraction peak and do not satisfy the 
phase function normalization condition. 


2.2.4 Polarization Characteristics 


Light beam is characterized not just by the intensity, direction of propagation and 
frequency. It also can be totally or partially polarized. The solar light incident on the 
top of terrestrial atmosphere is almost unpolarized. However, it can be polarized due 
to scattering processes in atmosphere and reflection/refraction by natural surfaces 
including snow. Polarization of incident unpolarized light by snow is rather weak in 
the visible. Any light beam can be characterized by the Stokes vector-parameter ር 
(I, О, U, V). The parameters J, Q, U, V describe the light intensity (T), degrees of 
linear (pe), circular (p<) and total (p) polarization 


Q? +U? V 
р = CELP ‚Р үр + pe. (2.125) 


I I 





The parameters (Q, U, V) also describe the orientation and ellipticity of 
polarization ellipse. The transformation of the incident Stokes vector parameter 


S (l, Qo, Uo, Vo) by a local snow volume can be presented in the following way: 
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— ANE 
5 = ТЕ 0, (2.126) 
where Ё is the scattering matrix, Y — (kr) 2, k = a, r is the distance to the 


observation point. The phase matrix P is proportional to the scattering matrix with 
the element Pi; equal to the phase function: 





; (2.127) 


where C, is the average scattering cross section. 
The normalized phase matrix elements are often used. They are defined as follows: 


pij = Pij/ Py. (2.128) 


In the case of randomly distributed nonspherical particles the phase matrix has 
the following general block-diagonal form: 


ыр е а рй (2.129) 
0 0 P33 Рд 
0 0 — P34 P44 


Therefore, Stokes parameters of scattered light can be presented in the following 
form: 


I = Piilo + Pi2Qo, 
Q = Pio + Р2 Оо, 
U = P33Uo + P34 Vo, 
V = —P34U9 + P44Vo. (2.130) 





It follows in the case of incident unpolarized light (/o = 1, Qo = Up = Vo = 0): 





I = Pi, Q = Pa, U = 0, V = 0. (2.131) 


This means that the scattered light becomes linearly polarized with the degree 
of polarization equal to the value of the normalized phase matrix element pio = 
Pis / Ру. Also one can easily establish physical meaning of other normalized phase 
matrix elements. Let us assume that the incident light is left handed circularly 
polarized (/6 = Vo = 1, Ug = Qo = 0). Then it follows from Eq. (2.130): 





I = Pu, Q = Pn, U = P34, V = Рд. (2.132) 
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The degree of circular polarization (DCP) is equal to pag = P44/ P11. Therefore, 
the element p44 is equal to the degree of circular polarization of singly scattered light 
under left handed circularly polarized light illumination conditions. Let us consider 
the illumination of a particle by the horizontally polarized light (/ = —Uo = 1, Vo = 
Qo = 0). Then it follows from Eq. (2.130): 





I = Pilo, Q = Рр, U = 0, V = Рм. (2.133) 


Therefore, the element p34 is equal to the degree of circular polarization of singly 
scattered light under horizontally polarized light illumination conditions. It shows the 
effectivity of linear to circular light polarization state transformation during single 
scattering events. In a similar way one can establish the physical meaning of other 
phase matrix elements. 

The normalized phase matrix elements of Voronoi aggregate and fractal particle 
models derived in the framework of geometrical optics approach are given in 
Fig. 2.13. The calculations have been performed at the wavelength 0.56 um, n = 1.31 
and the size parameter x — 311 (for Voronoi aggregates) and the characteristic length 
100 um (the side of initial tetrahedron) of a Koch fractal of second generation. It 
follows that the normalized phase matrix elements almost coincide for these diverse 
two models of ice grain shapes. One can see that such particles only weakly polarize 


normalized phasae matrix elements 





0 20 40 60 80 100 120 140 160 180 
scattering angle, degrees 
Fig. 2.13 Normalized phase matrix elements of ice crystals at the wavelength 0.56 ит (n = 1.31) 


calculated using fractal (tetrahedron of second generation with the side of 100 ит) and Voronoi 
(the size parameter 311) aggregate irregularly shaped particles 


2.2 Local Optical Properties of Snowpack 51 


incident unpolarized light (the degree of linear polarization рә is smaller than 10% 
with maximum polarization in the side scattering direction). The value of DCP p44 
decreases with the scattering angle. The value of the phase matrix element p34 
(linear to circular polarized light conversion probability) is small for all scattering 
angles. Such behavior is consistent with measured polarization characteristics of 
snow layers. 


2.2.5 Light Absorption by Polluted Snowpack 


Snow can contain various inclusions (soot, dust, algae, litter from the trees, etc.). 
Pollutants modify snow optical properties and appearance. It is often the case that 
the concentration of pollutants is too small to influence light scattering processes 
in snow. However, snow hardly absorbs light in the visible. Therefore, even small 
concentrations of pollutants influence local absorption coefficient of snow. Assuming 
external mixture of pollutants and ice grains, one may write for the snow absorption 
coefficient: 


kaps = kice + k?” (2.134) 


abs abs? 


ice 


where the ice absorption coefficient К, 


coefficient of pollutants P 


is given by Eq. (2.104) and the absorption 
can be presented as 


КРЛ = сур, (2.135) 


where c, is the volumetric concentration of pollutants and 


[2] 
[ሃዛ 





ур = (2.136) 


is the volumetric absorption coefficient (VAC) of pollutants. Here, (Уо) 15 the 


average volume of impurity particles and (с B | is their averag s section. For 


very small and also for very large ( as compared to the wavelength kly absorbing 


pol 


particles, the cross section (c abs 


| is proportional to the product of the bulk absorption 


coefficients of pollutants œ po; = xpa and (ол). Here, ጀ in the imaginary part 


of impurity complex refractive index m jo-7 pol — i Хро: Then it follows: 
Yp(A) = B pot pol (A). (2.137) 


In particular, it follows at the size of pollutants much smaller than the wavelength: 
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On poi 


(2 + mu) (2 Lnd Жз) + Хы 





Bol = (2.138) 


Assuming that the soot complex refractive index is 1.95—0.791 across solar 
spectrum, one derives: Вы = 1.3. 
One can derive the following simple equation for snow contaminated by soot: 


Карз = Cice Biceice + C pol B pot pol (2.139) 


lid in the visible and near infrared (at small values of the probability of photon 
ption). 


In general, the value of (с i | can be derived averaging the absorption cross 


section of pollutants C? ? with respect to their size and shape distribution. It follows 


abs 
“б Bs (a) f (a)da, (2.140) 


for the spherical particles: 
0 


where f(a) 15 the size distribution of impurities. One can see that the calculation 


of ic | А, requires information not only on the complex refractive index of impurity 
particles but also their size distribution must be known. This is often not the case. 
Therefore, various approximations of the VAC y, is used. In particular, one can 


assume that 


X -P 
Yp (А) = »eo(: ) . (2.141) 
Ао 
Опе can also define the mass absorption coefficient (MAC) of pollutants: 


My) = Z, (2.142) 


Here p is the density of an impurity particle (1.8 g/cm? for black carbon (BC) or 
soot, 2.5 glcm? for dust). The black carbon MAC M, (A) decreases from 6.8 r3 at 


0.4 um to 3.6 = at 1 um. It is equal to 0.7 me at 4 wm. 

The assumption of the external mixture of pollutants gives the lower boundary 
of the influence of pollutants on the snow absorption. As a matter of fact, some 
pollutants are contained inside snow grains (internal mixture). This leads to even 
greater influence of pollutants on snow absorption due to focusing properties of snow 
grains. The influence of pollutants on snow absorption is at maximum in the visible, 
where the maximum of the incident solar flux is located. It generally diminishes with 
the wavelength and can be neglected at channels, where ice strongly absorbs solar 
light (Dombrovsky and Kokhanovsky 2021). 
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Chapter 3 A) 
Radiative Transfer in Snowpack EM 


3.1 Radiative Transfer Characteristics 


The main task of this chapter is to derive simple analytical equations, which govern 
light intensity and fluxes inside snowpack and also at its boundaries. The spectral 
light intensity T, is defined as the amount of radiant energy dW in the wavelength 
interval [A, A + dÀ] transported across an element of area dS into a solid angle dQ 
oriented with an angle 9 relative to the normal to the area dS, within a time interval 
dt: 


h= се ብ. (3.1) 
cos 9#4О4Ааба:ї 
Therefore, if one derives the spatial and angular distribution of the light intensity 
I;, the amount of radiation energy dW transported through an element of area dS 
into a solid angle dQ in the spectral interval dÀ and time interval dt in snow can be 
calculated as 


dW = І, cosódQdAdSdt. (3.2) 


It is useful to introduce irradiance as the amount of energy transported through 
an element of area dS in in the wavelength interval dÀ and time interval dt 





dW 

= . (3.3) 
dAdSdt 
The radiant power 

aw 

аф = — (3.4) 
dt 
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is the amount of transported energy in unit time. It is measured in watts (J/s). The 


irradiance is measured in [w que nm !] а 
An important point is that ће monochromatic irradiance an be calculated 
as [see Eqs. (3.2), (3.3)]: 


Ka j= (s. 9) 608942, (3.5) 


where the light intensity T, (7, 9) is measured in [w m^? . 8171. am !]. Itis a 


Я Я > ‚ К dol E A 
function of coordinates r, propagation direction (2 and time f: 


h=(, El (3.6) 


where r is the position vector incorporating coordinates, g is the unit vector in 
the direction of light propagation, and f is the-time. The intensity does not change 
with time assuming that the properties of the light propagation channel (say, a turbid 
medium) and light source do not change in time (stationary media and light sources). 
One should account for the time dependence of light intensity, e.g., in the case of 
problems related to the illumination of snow by a pulsed (say, laser) light source. In 
this book we shall consider only stationary light fields and turbid media. 

For the medium stratified in plane-parallel horizontally homogeneous layers such 
as a snowpack, the intensity depends on the vertical coordinate z and the propagation 


direction g : 
— 
h=h(z 0) = AG o. 7) 


Here the direction g = {u, ф} is defined in spherical coordinates through the 
cosine u of the zenith angle 0 and azimuthal angle o, и > 0 for downwelling 
radiation and u < 0 for the upwelling radiation. Note that 

dQ = dudg. (3.8) 

If the light intensity is independent on direction, then it is called isotropic, that is 

L. = L.G2). (3.9) 
If the intensity is integrated over all wavelengths, it is called the integrated intensity 


and given by 


oo 
r= | hai. (3.10) 
0 
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The irradiances for the downward and upward propagating light are defined as 


2л 1 
ri f | ne ሥሪ (3.11) 
0 0 
and 
2x 0 
rto = | | һи, ийиде. (3.12) 


0 ፦1 


Quantities (3.11), (3.12) are called downwelling and upwelling light fluxes, 
respectively. The integrated irradiance over all wavelengths reads as 


oo 
P= | ња. (3.13) 
0 


Yet another quantity of interest is the radiation energy density. It is defined as 


1 
оу = | nag, 3.14) 
C 
4л 


where c is ће light speed. The actinic flux 


ቁዋ, = co; (3.15) 


is needed for the calculation of absorbed radiation power at a given location in snow. It 
particular, this radiative characteristic can be used for the studies of photodissociation 
rates of various molecules in snow. 

One can see that the light intensity T, is the main radiative transfer characteristic, 
which defines all other optical properties (e.g., radiant power, irradiances, actinic 
fluxes, radiation energy density, etc.). In the next section we shall introduce the 
equation, which can be used to calculate the light intensity for plane parallel hori- 
zontally homogeneous snowpacks uniformly illuminated at a top by a plane-parallel 
light beam. 

It should be pointed out that the intensity 7; is related to the photon distribution 


function EP by a simple equation: 
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I = chv f, (r, 9), (3.16) 


where є = hv is the energy of a photon jn) the Planck constant and v = c/A is the 
E > P 
frequency. The integral N(7) = f fi(r, $2)d €? gives the number of photons at the 
4л 


position specified by ће vector 7. Equation (3.16) follows from the comparison of 
Eq. (3.2) with yet another representation of the radiant energy: 


dW = edN, (3.17) 


where 


dN = ^(f. 8) cos 9dQdAdV (3.18) 


is the number of photons with wavelengths in the interval [A, А + dA] located in time 
t in a volume dV = cdtdS near point r and moving within a solid angle dQ around 


direction б. It follows from Eqs. (3.14), (3.16) for the radiation energy density: 
Pus m | fl? 9)ао. (3.19) 
4л 
Therefore, the ratio 


Pr 
= = 3.20 

hv (3.20) 
gives the number of photons with the frequency v in unit volume of a snowpack. 
This parameter is of importance for the snow photochemistry. 


3.2 Radiative Transfer Equation 


The snow radiative characteristics can be derived using the radiative transfer equation 
(RTE) for the diffuse light intensity 7. This equation can be written in the following 
form for the plane-parallel snow layers (Chandrasekhar 1950): 


2л л 
dI 0, x ያ г a 1 
oso OO (о, e) e fag fav p(9'. o! > 9,9) Q^.) 
dt 4л 
0 0 


+ 2 p(o, ф > 9, 9) Foe >. (3.21) 
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Here wp is the single scattering albedo, т is the optical thickness, p (9', ф' > 9, ሀ) 
is the phase function, ፲ Fo is the incident solar flux at the unit area perpendicular 
to the light beam, #0 is the solar ineidenee angle, фо is the solar azimuthal angle, 
9 is the viewing zenith angle, and ф is the viewing azimuthal angle. This integro— 
differential equation enables the calculation of the diffuse light field intensity at any 
level and direction of propagation inside of the snow layer and also at the upper and 
lower boundaries of a snowpack. The calculation must be done taking into account 
boundary conditions. The simplest approximation is to assume that there is no diffuse 
light incident on the snow layer from above and also from below (underlying black 
surface). Such a condition is close to reality, e.g., in Antarctica, where atmospheric 
optical thickness is very small and diffuse light coming from the sky can be often 
neglected (say, at near infrared region of electromagnetic spectrum). For polluted 
atmosphere, the account for the diffuse light incident from the sky is of importance. 
In this case one needs to specify the vertical distribution of local optical properties of 
atmosphere—snow system (Aoki et al. 1999) and solve Eq. (3.1) numerically (say, 
in the framework of the discrete ordinates technique). 

The direct light beam intensity unlike diffuse light intensity does not require 
complex numerical calculations. Namely, it follows: 


lie = Fo ew(- s oct — Оу), (3.22) 
0 


> > > > 
where the delta function ፅ( (2 — © ዐን signifies the propagation direction ( 62 = Qo 
of the direct beam. It follows that the direct beam attenuates faster for larger solar 
zenith angles #0. 

The integral in Eq. (3.21) accounts for multiple scattering effects. The single 
scattering contribution to the reflected and transmitted light can be derived from 
Eq. (3.21) ignoring the integral term. The remaining differential equation can be 
solved analytically. The answer is: 


x 


Fo —X / / T 
I(x) = 1° | w(x )p(x )exp(- (q — Dx), (3.23) 
y 
where x — ፲/ 665 9,4 —cos 9 / cos до. The constant V can be found from the 


boundary conditions. Let us assume that the diffuse light intensity incident both 
on the top and bottom of a scattering layer is equal to zero: 


LG = 0) = 0, соѕд > 0, ha = хо) = 0, cos? < 0. (3.24) 
Here arrows show the direction of light propagation, x = 0 corresponds to the 


upper boundary, x — xo corresponds to the lower boundary of the scattering layer 
with the optical thickness zo and xo = To / cos 9. Let us assume that the local optical 
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properties do not change with the depth. Then the integral (3.23) can be evaluated 
analytically: 


eap ду —x _ gt]. 


(х) = oc ፦ [e= 0D _ 4х}, (3.26) 


where 0 is the scattering angle and we have accounted for the boundary conditions. 
It follows at the boundaries: 





0) F | 
1, (хо) = ne “ከሁ” eo), (3.27) 
1,0) = > T [e 79» — 1]. (3.28) 


Let us introduce reflection and transmission functions: 


mh (0,@,z = 0) 

















R , 3 , = 3.29 
(и, Ф, Ko, Фо) Fok (3.29) 
л1,(Ф#„ф,т=т 
ТУК с = (3.30) 
Foé 
where £ = cos 9o. Then it follows: 
— ማዘሆ፥ህ)”] 
R= 1 — ех + Tol? (3.31) 
4( +m) He a) 
БЕ ji 
T= exp exp , (3.32) 
4@ — n) n š 
where 7 = |cos 9| and cos0 = —£n + y (1 — £2)(1 — n?) cos(q — qo). 
Equation (3.32) has the following form at £ = n: 
@oto p (0) то 
T= ДЕ? | Е | (3.33) 


The contribution of single scattering events to the transmission can be neglected 
for very large values of ro typical for snow layers. However, the contribution of singly 
scattered photons to the reflection function can not be neglected even for the case 
of semi-infinite layers and this contribution is given by the following equation [see 
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Eq. (3.31)]: 


_ wp) 


= ሚማ 3.34 
4(@ + m) om) 


SS 


This equation shows that not only integral scattering and absorption characteristics 
but also snow phase function must be accurately modelled to predict snow reflective 
properties (e.g., angular dependence of snow reflection function). The contribution 
of multiply scattered photons to the reflection function is discussed in Sect. 3.4. 


3.3 Light Field in Deep Layers of Semi-Infinite Weakly 
Absorbing Snowpack 


Let us assume that the medium under consideration is a semi-infinite one and the 
light intensity at a large optical depth inside the medium is to be derived. In this case 
one may neglect the contribution of the direct light beam. Then it follows: 


1 


210 = —I(u) + > J I(u', u)durp(cos 0), (3.35) 


፦1 


where J (jz) is the diffused light intensity in the direction specified by the angle 
0 = acos(u) from the normal to the layer and 


2л 


7(605 ፀ) = x= “መ (3.36) 
2л 
0 


is the azimuthally averaged phase function. 
The intensity in deep snow layers does not change its angular pattern and is 
attenuated according to the exponential law (Sobolev 1975): 


I = Ci (u) exp(—kt), (3.37) 


where k is the diffusion exponent, C is the constant. 
Let us substitute this solution to Eq. (3.35). Then it follows under assumption of 
isotropic scattering (p — 1): 


1 
1—Ки 





i(u) = : (3.38) 
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where we used the normalization condition: 
= Sitaw e (3.39) 


Substituting Eq. (3.38) to equation Eq. (3.39) given above makes it possible to 
derive the following equation for the diffusion exponent: 


wos(k) = 1, (3.40) 
where 
1 1+К 
К) = — In ——. 3.41 
MAS TEE Ga 


Let us take into account that it follows for weakly absorbing media (k — 0): 


5 КЎ 
Ind+k)=k-—+4+-—. (3.42) 
2 3 
Therefore, one derives: 
k2 
s(k)=1+ 3 (3.43) 
and [see "Loos 
1+ E NM. (3.44) 
3 1- | 
where В = 1 — o is the probability of phot sorption. Finally, we derive for 


weakly absorbing media (6 — 0) from Eq. (3 (at p = 1): 


к = 3g. (3.45) 


This equation shows how the diffusion exponent depends on PPA for the case of 
weakly absorbing isotropically scattering media. Therefore, one derives at p — 1 and 


k > E] 
и) = 1 + /38и. (3.46) 


The phase function of snow is not isotropic and asymmetry parameter is close to 
0.75 in the visible. Therefore, it is of importance to consider more complex phase 
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functions. This can be done solving Eq. (3.35) numerically for the phase func- 
tions specified in the previous Chapter. However, our task is to derive approximate 
analytical equations. Therefore, we consider now the case of the following phase 


function: 


p(0) = 1 + xi cos 0, (3.47) 


where x, = 3g, g is the average cosine of the scattering angle. Then one derives 
assuming that the value of / is given by Eq. (3.37): 


z 00 i Nefa! 7 
a -kmi = > I (1 хин? u!) dut (3.48) 
BET 


where we accounted for the fact that the azimuthally averaged cosine of scattering 
angle can be presented as (cos Ө)„ = ши’. Therefore, we derive: 








s b(u) 
i) = D (3.49) 
where 
1 
00 IS. Й 
buo = 2 | ቦ+።ሠዞባሠ (3.50) 
E 
or 
"(+ib (u) 
x 
Ь(и) = >f — Р dy. (3.51) 
21 


It is clear from Eq. (3.50) that b(u) can be presented in the following form: 





b(u) = bo + ри, (3.52) 
where 
^b ^d / i ЙЛ ) 
Бу = = | SHOE Ж а መ [t (и аш. (3.53) 
2 ]—kyw 2 ]—kyw 


—1 =} 
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It follows from Eqs. (3.52), I 


bo + biu = bo jo + biji + box jiu ተ bixi jau (3.54) 
or 
bo = bojo + biji, bi = boxi ji + bixi j, (3.55) 
where 
i / nq / 
Јһ = Фо ("du (3.56) 
2] 1-kw 


or after integration 


А . e(s—-l), j 
jo = os, ji = — Т pet) (3.57) 





Therefore, we have the following equation for the determination of both bo and 
bi: 


ab = 0, (3.58) 
where the elements of the matrix à can be represented in the following way [see E 
(3:35).-G.37)]: a 


1 1—5 1 s—1 (3.59) 
=——-S, a =a = , a መመመ — . . 
ai m $ 12 21 k 22 UH B 








This equation has a solution, if the determinant 


D = апа — аро (3.60) 


is equal to zero. Therefore, we have: 


] — a 
k2 





1— Wo 
w| 1+ xi P S — WX] = 1. (3.61) 


Equation (3.61) can be used to determine the diffusion exponent at any values 
of the pair (x1, ወ0)› Let us assume that the medium under consideration is weakly 


absorbing. Then the probability of photon absorption В = 1 — wọ 0, 5 > 1+ E 
and, therefore, one derives: 
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k= y 3B — ғ), (3.62) 


which is reduced to the value derived for the case of isotropic scattering in the case g 
= (as it should be. Importantly, although this expression for the diffusion exponent in 
weakly absorbing medium has been derived for the case of a specific phase function, 
it can be used for any phase functions under assumption that 6 — 0, which is valid 
in particular for snow in the visible and near infrared. Assuming that g — 3/4 for 
snow in the visible and near infrared, we derive: k = 4/38 /2 for snow in the visible 
and near-infrared. One can see that the diffusion exponent is two times smaller for 
snow as compared to the case of isotropic scattering under approximation studied. 
Therefore, the radiation penetration in snow is larger as compared to the case of 
isotropic scattering for a given value of probability of photon absorption. 
It follows for the function i (u) 


1+é&u 
1—kp’ (3.63) 





(ሠ) = bo 


where £ = b,/bo. This function can be multiplied by any constant and still be a 
solution of Eq. (3.35). Let us assume that bọ = 1 and £ = bi. Then it follows from 
Eq. (3.55): 





ደ = 317. ( + i) (3.64) 
or 
£= = (3.65) 
~t— Te i 
where П = (s — 1)wọxık™!. Let us assume that В — 0. Then one derives: 
k 
a (3.66) 
1 — goo 
and 
i(u) = 1 + Зѕџ, (3.67) 
where 
5 = В (3.68) 
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is the similarity parameter. One can see that the ratio of the probability of photon 
absorption to the diffusion exponent is a key parameter of the radiative transfer in 


snow. Also we can write: 
k 
5 = f = д (3.69) 
3⁄1—g) 23(1—9) 


The parameters (k, s) can be used to derive the intensity of light in deep layers 
of weakly absorbing media with the phase function given by Eq. (3- Namely, it 


follows: 
I(u) TI 351) exp(—kT) (3.70) 


as В — 0. This equation can be also written in the following form: 





Iw) TI 3s) exp(—yh), (3.71) 
where 


y = kkext (3.72) 


is the so-called asymptotic flux attenuation coefficient (AFEC) and h is the geomet- 
rical depth. It follows from Eq. (3.71) that the ratio ¢ of light intensities/fluxes at two 
depths (hi, h>) in deep layers of snow can be presented as 


t = exp(—y A), (3.73) 


where A = h> — hi. Therefore, the AFEC can be easily measured experimentally: 


_ I(F(hi)/(F(ho)) 





(3.74) 


Here, F (h j) is the flux at the depth h;. Equation (3.74) is also valid for light 
intensities in a given direction as well. The measurements of light intensities at 
two directions at a given depth in snow can be used to estimate the parameter s 
[see Eq. (3.71)]. Although Eqs. (3.62), (3.70) have been derived for a particular 
phase function, they play a fundamental role in the radiative transfer in weakly 
absorbing media with arbitrary phase functions. As a matter of fact, they are valid 
approximations for the snow in the visible and near infrared. Our derivations show 
that the most important single local optical characteristic of snow is the similarity 
parameter s. This parameter can be used to find the diffusion exponent k = 3(1 — 
g)s (under assumption that g is known). Also light intensity distribution in deep 
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snow layers [see Eq. (3.70)] can be derived. The parameter s depends on the ratio 
of absorption laps = 1/kabs and transport extinction ያሉ = 1/(1 — g)k¿a (valid at 
kabs << Куса) lengths: 


ltr 
3labs | 





(3.75) 


= 


Also it follows: /,, / labs = 3s?. Generally speaking, the two local optical parameters 
(y, s) are needed to describe the light propagation in snow layers. The value of the 
diffusion exponent can be derived from the AFEC in case the snow extinction coef- 
ficient kext is known [see Eq. (3.72)]. We shall see later that the similarity parameter 
s can be derived not only from the analysis of the angular distribution of transmitted 
light in deep snow layers but also it can be obtained from the measurements of snow 
diffuse reflection coefficient under diffuse illumination conditions. This means that 
the transport extinction length can be also derived: 


3 
per (3.76) 
y 


In addition, the absorption length is given by the following equation: 


A 
ማቸ a (3.77) 


This opens ways to the reflectance spectroscopy of snow (in particular, the 
determination of pollution type and load from measurements of snow reflectance 
spectra). 


3.4 Reflection of Light from a Semi-Infinite Snow Layer 


3.4.1 Nonlinear Integral Equation for the Reflection 
Function 


Baud > can be used for any optical thickness. In many cases snow can 
be considered as a semi-infinite medium. Then it is useful to use the following 
representation for the reflection function: 


R= Rss + Rms- (3.78) 


The value of Rss gives the contribution of single scattering [see H... and 
Rms describes the contribution of multiple scattering to the reflected light intensity. 
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The value of Rms can be calculated using Eq. (3.1) or alternatively using this nonlinear 
integral equation valid for semi-infinite turbid media (Ambartsumian 1943): 


ex (š L() + nL(5) + £&nM E, m) 





RoR, (3.79a) 
£ +n 
where 
1 1 2x 
LG) = — | / p(n, бай ФУКО, Ф', E, edi dg, (3.79b) 
0 0 
| 1 2x 
56) = 2. | | ne wa CORO On dared’, (3.790) 
0 0 
2л 1 
1 ያ là { ያ 
М(8) = 25 | 4 ат R(n', Ф', n. Ф) 
0 0 
2л 1 
х faw [ovo е. n”, ф")К(п", Ф", E, фо). (3.794) 
0 0 


Опе can see that L and M depend on the reflection function А. There are different 
methods to solve the nonlinear integral equation described above. One possibility 
is to use the successive order approach, when Rms is estimated from Eq. (3.792) 
assuming that R = Rss in the first approximation. This enables the calculation of 
the double scattering approximation for the reflection function, which can be used to 
derive the triple scattering approximation, etc. Such series are converging quite fast 
at small values of wọ. Otherwise, a lot of terms are required to reach a convergence. 
The functions L and M depend just on the phase function and the single scattering 
albedo and can be parameterized using numerical radiative transfer calculations of 
the reflection function. 

Let us assume that isotropic scattering occurs in a nonabsorbing semi-infinite 
medium. Then both plane r, and spherical r, albedos are equal to 1.0. They are 


defined as 
1. 1 
r (Š) = 27 RE, n)ndyn, r; = 2 ወ (3.80) 


where 


ze 1 2x 
R(E, n) = > / RG, n, #)ሪሃ- (3.81) 


3.4 Reflection of Light from a Semi-Infinite Snow Layer 69 


It should be pointed out that the reflectance R is defined as the ratio of intensity of 
reflected light from a given target to the reflectance of absolutely white Lambertian 
surface. The isotropic scattering at single scattering albedo equal 1.0 (no absorption) 
is very close to the case of a Lambertian white surface. Therefore, one can assume 
that R is close to one. Then one can use the following approximation: 


one 


RE, dg = "E (3.82) 


This means that L — 1/2, M — 1 [see Eqs. (3.79)] and it follows for the isotropic 
scattering from Eq. (3.79) at the single scattering albedo wo = 1: 





+n+2 
M (3.83) 
2( +n) 
Therefore, one derives for the case of isotropic scattering at oo = 1: 
1 1 
patpi (3.84) 


252508. AEG» 


We assume that the dependence of L and M on the angles can be neglected in the 
first approximation also for other types of turbid media. Then it follows: 


E 11(፪ +n) + 13፪# + 173 
di 2(§ +n) 





(3.85) 


and the constants Г; (j = 1, 2, 3) can be derived from fitting te this equation to 
the radiative transfer calculations. It follows from the discussion given above that 
Г! ጊጋ = 3 = 0 for isotropic scattering. It follows (Kokhanovsky 2005) for 
the al о ይክ phase function: Г, = 0.593, Г = 2.579, and Гз = 0.624. 
Therefore, the following approximation for the nonabsorbing snow reflectance can 
be used: 


. ME +n) + Г»ёт + Ts + p(0)/2 


Ro = 7 3.86 
0 XE +m (3.86) 





where p(0) is the snow phase function. One can use the parameterizations of the 
phase function as given by Eqs. (2.114) and (2.115). The reflection function of the 
semi-infinite turbid layer calculated using RTE (solid lines) and Eq. (3.86) (dashed 
line 1) for turbid media composed of nonabsorbing spherical water droplets (water 
cloud) and ice fractal particles (snow) as the function of the incidence angle at 
nadir observation and the wavelength 550 nm is shown in Fig. 3.1. The spherical 
water droplets are characterized by the effective droplet radius 6 microns and the 
coefficient of variance of gamma size distribution 38% (Deirmendjian 1969). The 
phase function of the fractal particles have been calculated using the fractal model of 
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Fig. 3.1 The reflection 1.5 - - - - - - - 
function of the semi-infinite ኮ 
turbid layer calculated using 
RTE (solid lines) with the 
code developed by 
Mishchenko et al. (1999) and 
Eq. (3.86) (dashed line 1) for 
turbid media composed of 
nonabsorbing spherical 
water droplets (water cloud) 
and ice fractal particles 
(snow) as the function of the 
incidence angle at nadir 
Observation and the 
wavelength 550 nm. The 
refractive index of water and 
ice was taken from Hale and 
Querry (1973) and Warren 
(1984), respectively. The 
dashed line 2 corresponds to 
Eq. (3.86) with р(Ө) = 0 
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ice grains (Macke and Tzschihholz 1992; Macke et al. 1996). It follows that Eq. (3.86) 
can be used to describe the angular behaviour of the snow reflection function at 
channels, where light absorption effects by snow can be neglected. One can see that 
the reflection function for the case of fractal snow grains is smooth and closer to the 
Lambertian case (R — 1) as compared to the water cloud model. Also there is no 
peak for the reflectance of turbid media with fractal particles in the rainbow region 
characteristic for water droplets. Therefore, one concludes that the assumption of 
spherical ice particles can lead to substantial errors in the interpretation of angular 
light intensity dependence for the case of light reflected from snow layers. This is not 
so crucial for the light transmitted by thick snow layers because of randomization of 
photon directions. 


3.4.2 Reflection Function of Semi-Infinite Weakly Absorbing 
Snow Layers 


Letus consider the case of weakly absorbing semi - infinite snowpack. The reflectance 
can be presented as (Rosenberg 1962; Kokhanovsky 2002): 
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oo 


RB) = 33a. — ይ". (3.87) 


n=1 


where В = 1 — оо is the probability of photon absorption by unit volume of snow. 
In the case of nonabsorbing snow it follows: 


oo 
RO) = Ma. (3.88) 
n=1 
Let us define the ratio: 
9t = RCB)/ R(O). (3.89) 
Then we have: 
oo 
> ax, (1 m В)" 
R= ።፤ (3.90) 


ጩ 
> a, 


n=1 
Expanding the nominator of Eq. = we derive at small B: 


2 3 
ሸ51-ቆ<ክ>+9. < > -Ê < m > +... s< exp fn) >, (3.91) 


where 


ап 


У fun, (exp(—Bn)) = » /,ехр(-Вп), f, = = (3.92) 
n=1 n=1 Э аһ 


п=1 





(t) 


and we assumed that that n(n — 1) ~ n?, n(n — 1)(n — 2) ~ ዘ”... іп our derivations. 
This is possible because P is close to zero and the number of scattering events is high. 
For the same reason we have: 


= x —Bn)d 3.93 
< exp(—fn) > J f(n) exp(—pn) "E] (3.93) 


This integral can be evaluated assuming the function f (n). In particular, it follows 
from the random walk theory that that the probability of a particle (photon) appearing 
ata given place, time, and direction after large number of iterations n can be presented 
as (Chandrasekhar 1943): 
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fO) = (oer exp(—v/n). (3.94) 


The parameter u depends on the process studied. The substitution of Eq. (3.94) 
into Eq. (3.93) gives: 


9t = exp( 2 / uf). (3.95) 
Therefore, we can write: 
R(B) = Roexp(—y ОВ), (3.96) 


where Ry = R(0), Q = 4 v. This equation has been proposed for the first time 
by Rosenberg (1962). It shows that the spectral snow reflectance depends on the 
square root of the probability of photon absorption 6. The paramete pends on 
the scattering and not on absorption processes and, therefore, one may assume that 
it does not depend on the wavelength for snow composed of large snow grains in 
contact. It should be pointed out that we have derived Eq. (3.96) not using the radiative 
transfer equation. Actually, the applicability of RTE for closely—packed media such 
as snow is in question. It should be pointed out that Eq. (3,79) is very general and 
can be applied to many geophysical media such as snow, leaves, etc. The main 
shortcoming of this equation as applied to natural snow is the fact that it has larger 
errors in areas close to forward scattering (glint) and backward scattering, where 
snow has enhanced brightness, e.g. due to oriented snow crystal surfaces or thin ice 
films on the top of snowpack. This equation can be applied only for the horizontal 
smooth snow surfaces. In particular, a modification is needed for a sloppy terrain 
(Picard et al. 2020). Zuravleva and Kokhanovsky (2011) has studied the influence of 
surface roughness on the radiation reflected from snow surfaces. It was found that 
the albedo of a snow layer is reduced (in particular, in the infrared region), if 3D 
effects are taken into account. 

The value of Q can be related to the asymmetry parameter g of ice grains 


using asymptotic results of RTE y4 at small values of the probability of photon 
absorption. Then it follows from “733: 
R(B) = Ro(1 — y ОВ) (3.97) 


and also from analytical radiative transfer theory (Sobolev 1975): 


R(B) = Ro — yu (uou (u), (3.98) 


where 
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y=4 i-a (3.99) 


It follows that y — 4s (see Table 2.1). The escape function can be found from the 
following equation (Sobolev 1975; Minin 1988): 


3 
uo (i) = 4 Ho + #(ሠ0)))› (3.100) 
where 
t u 1 1 
Vuo) = 2 | Rolo, шиаи, Ro = 7— | Коби, но, Фф аф. (3.101) 
0 0 


For weakly absorbing semi-infinite media one may expect that Ro is close to one 
and the dependence of W on the angle can be neglected [see Eq. (3.100)]. Therefore, 
the escape function can be approximated by the linear function of uo. 

The following normalization condition holds (Sobolev 1975): 


1 
2 f uo (o) Mod Ho =1. (3.102) 


Also it follows for a special case of isotropic scattering (Sobolev 1975) 


1 
ያ uo(Ho)d Ho = 2uo(0). (3.103) 
Let us assume that the escape function can be represented as 


uo(uo) = ።0(0)(1 + cuo). (3.104) 


One derives after substitution of Eq. (3.104) in Eq. (3.102), (3.103) that c = 2 and 
ሀ0(0) = 3/7 for the case of isotropic scattering. This means that 


3 
uo(uo) = 70 + 270) (3.105) 


for the case of isotropic scattering. The escape function ዘ(/0) has been tabulated 
by Yanovitskij (1997) for the Henyey-Greenstein phase function with the asymmetry 
parameter g — 0.75 close to that of snow. The results are given in Table 3.1. The data 
shown in Table 3.1 can be approximated by the following simple equation: 


3 1+ ኣ/# 
toto) = Sito —— (3.106) 
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Table 3.1 The escape function for the Henyey-Greenstein phase function at g — 0.75 (Yanovistskij 
1997) 





o |0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 





ио |0.333 |0.484 0.95 | 0.692 |0.783 | 0.869 | 0.952 | 1.033 | 1.113 |1192 | 1.27 



































Both Eqs. (3.105) and (3.106) satisfy the normalization condition (3.102). They 
have a similar accuracy at o > 0.2. However, Eq. (3.106) gives better results at 
smaller values of ио (see Fig. 3.2). 

Comparing Eqs. (3.97) and (3.98), one derives: 


16 f? 
Q= ae (3.107) 
30 — 8) 
where 
f= и(ио)и(и) (3.108) 
Ко 
Therefore, ме have finally (Хесе et al. 1991): 
Roo = Ro% ехр(—и(и)и(ио)у/ Коо), (3.109) 


where we introduced the symbol оо to underline that we consider a semi-infinite 
turbid layer and y is given by Eq. (3.99). Essentially, Eq. (3.109) reduces the problem 
of the calculation of the reflection function of weakly absorbing turbid layer to that 
of a nonabsorbing layer (Ros) with the same phase function. The reflection function 


Fig. 3.2 The dependence of 
the escape function on the 
cosine of the solar zenith 
angle for the 
Henyey-Greenstein phase 
function at g = 0.75 
calculated using exact 
radiative transfer theory 
(symbols), Eq. (3.105) 5 
(dashed line) and Eq. (3.106) 
(solid line) 
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Fig. 3.3 The dependence of the reflectance of a semi-infinite snowpack on the cosine of the 
solar zenith angle at various values of single scattering albedo (SSA) and nadir observation. The 
Henyey-Greenstein phase at g = 0.75 has been used (exact results: crosses, approximation: boxes) 
(Kokhanovsky et al. 2019) 


Ros; depends only on the snow phase function and observation geometry. It can be 
easily parameterized as discussed in the previous section (see also Fig. 3.1). 

The accuracy of Eq. (3.109) is demonstrated in Fig. 3.3 as compared with exact 
solution of RTE at the Henyey-Greenstein phase at g — 0.75 and various values of 
single scattering albedo (SSA). One can see that the analytical approximation is very 
accurate and can be used for the solution of many applied problems. The value of Ros; 
has been calculated using the numerical solution of the radiative transfer equation. 
The error of the approximation is increased if an approximation (instead of exact 
value) for the function Ros; is used. One can see that reflectance is smaller at the Sun 
located at the horizon (small values of cosine of incidence angle) as compared to the 
case of the illumination of snowpack along nadir direction (at the nadir observation). 
The reflectance of semi-infinite snow layers is in the range 0.2-0.3 at the single 
scattering albedo 0.95, which is quite low. This means that the probability of photon 
absorption larger than 0.05 leads to very small values of reflectance. Then the weak 
absorption approximation is not valid and snow is not a bright object anymore. The 
reflectance decreases with the solar zenith angle for nonabsorbing snow and also 
for a polluted snow with small values of p (see Fig. 3.3). For larger values of p, the 
dependence of reflectance on the solar zenith angle becomes nonmonotonic function 
with a maximum located around SZA — 60 degrees. This is due to the fact that photons 
have more chances to be absorbed in the medium in the case of nadir illumination. 
Therefore, the tendency of high light reflectance in the nadir observation direction 
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characteristic for nonabsorbing media (as compared to the case of Sun at the horizon) 
is less pronounced or not valid anymore in the case larger values of probability of 
photon absorption (see Fig. 3.3). 


3.4.3 Snow Albedo 


Equation (3.109) can be used to find the plane albedo. For this, we shall use weak 
absorption limit of approximation (3.109): 


Roo = Roo; — и(и)и(ио)у (3.110) 


Then it follows from Eq. (3.80), oka (3.110): 


Гроо = 1 — и(ио)у, (3.111) 


where we accounted for the fact that r pœ = 1 for nonabsorbing media by definition. 
We shall use the exponential approximation similar to Eq. (3.109) to extend the 
applicability of the derived approximation with respect to the value of the parameter 


y: 
l'poo = exp(—u(4o) y). (3.112) 


Using similar arguments we derive the following approximation for the snow 
spherical albedo: 


Гуоо = exp(— у). (3.113) 
Equation (3.113) can be also derived using Feynman path integrals technique 


(Perelman et al. 1994; Feynman and Hibbs 1965). In particular, it follows (Perelman 
et al. 1994): 


n 
зоо = exp} —2 EX. (3.114) 





where Y is the constant, which does not depend on wy. Equations (3,110) and (3144) g 
coincide at wọ — 1 an = 4/3. 

It follows from e т TAN (3.410) (Zege et al. 1991; Kokhanovsky and Zege 
2004; Kokhanovsky et al. 2005a, 2007, 2018): 


Roo = Ror. 


soo? 


(3.115) 
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Рроо = Г”. (3.116) 
It follows from Eq. (3.116): 
F'poo (401) = rroo (402). (3.117) 


where к = ።(//01)/#(/፤02)› Moi is the cosine of the solar zenith angle до, zoo is the 
cosine of the solar zenith angle доз. This equation makes it possible to determine the 
snow plane albedo at any solar zenith angle, if the plane albedo is known (measured) 
at only one solar zenith angle (under the assumption that snow microstructure is not 
changing during the day). 


Equation (3.115) can be used for the determination of spherical albedo from the snow 
reflection function measurements: 


Fso = Uu Rue). (3.118) 
Then the plane albedo can be derived fro q. (3.116). The accuracy of 
approximations for the plane and spherical albed shown in Figs. 3.4 and 3.5. 
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Fig.3.4 The plane albedo as the function of the cosine of the solar zenith angle at various values of 
single scattering albedo for the semi-infinite snow layer with the asymmps=eparameter g = 0.75. The 
exact results are present ether with two approximations [boxes-H 2492), circles-calculation 
of plane albedo using s» (Kokhanovsky et al. (2020) ም 
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Fig. 3.5 The spherical albedo derived using Eq. (3.90) and exact radiative transfer calculations as 
function of the parameter y — 4s (Kokhanovsky et al. 2020) 


The spherical albedo у; can be used to determine the semi-infinite snow absorp- 
tance under diffuse illumination conditions: П = 1 — Fsoo. It follows in the limit of 
weak absorption: IT = y. This makes it possible to establish the physical meaning of 
the parameter y. Namely, this parameter provides the value of snow absorptance at 
В > 0. 

The parameter y can be related to the size of ісе grains in snow. In particular, it 
follows as 8 — 0 [see Eq. (2.109)]: 


у = Моб, (3.119) 
where ደ = edep is the effective ie» m length (EAL), which is proportional to 
the snow grain effective diamet d 

аа (3.120) 
pe ===, . 
9t = 8) 


The parameter ¢ depends on the shape of particles. However, its dependence on the 
wavelength and the size of particles in the spectral range 0.4—1 jum can be neglected 
(see the results presented in the previous Chapter). Therefore, we conclude that the 
EAL can be derived from spectral snow reflectance measurements even, if the shape 
of snow crystals is not known in advance. The determination of the effective grain 
diameter requires the parameter £, which depends on the shape of ice crystals/snow 
type. 
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It follows from Eqs. (3.113), (3.119): 
soc (À) = exp /a (0), (3.121) 


where a (A) = SAO) is the bulk ice absorption coefficient given in Appendix. 


Because the function о (À) is well known, we can conclude that the spectral diffuse 
reflection (under diffuse illumination conditions) of a flat clean snow layer depends 
on just one parameter, namely-the effective absorption length L, which can be derived 
from the measurements of the spherical albedo at a single wavelength in the near 
infrared (say, around 1 jum). The same is true for the plane albedo [see Eq. (3.112)]. 
The accuracy of this approximation is demonstrated in Fig. 3.6. It follows that the 
proposed approximation can be used in the spectral range 300-1200 nm with a high 
accuracy. The accuracy decreases at larger wavelengths due to the fact that the error 
of the weak absorption approximation increases and also due to the fact that one 
must account for the different penetration depths of radiation to different layers of 
vertically inhomogeneous snow. Equation (3.121) does not take into account neither 
surface roughness no vertical snow inhomogeneity effects. The slight discrepancy 
of measurements and theory around 400 nm can occur due to the fact that our theo- 
retical model ignores possible contamination of snow by various impurities. The 
uncertainty in the imaginary part of the complex ice refractive index (see Appendix) 
can play a role as well. 


* measurements 
theory 


spherical albedo 
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Fig. 3.6 The intercomparison of calculations using Eq. (3.121) (line) with experimental measure- 


ments (symbols) performed in Antarctica (Grenfel et al. 1994). The experimental results are derived 
using averaging of data (19 spectra) for two stations (Vostok and South Pole) 
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3.4.4 Snow Broadband Albedo 


In some cases the broadband albedo (BBA) is of interest. It is defined as 


À2 
/ “ዕ)ፆዕ.)ዐ› 


At 


ፈን ጮኢ (3.122) 
f FO)dX 


ài 


where F (À) is the solar flux incident on the snow surface and r (À) is either spherical 
or plane snow albedo depending on the illumination conditions. Clearly, for a cloudy 
sky, г (A) in Eq. (3.122) is the spherical albedo. For a cloudless clean atmosphere, r (À) 
is the plane albedo. We show the spectral solar flux F (A) for cloudless atmosphere 
at various solar zenith anlges (SZA) in Fig. 3.7. The results have been obtained 
using the SBDART (Ricchiazi et al. 1998) radiative transfer model. The absorption 
features of various atmospheric gases are clearly seen in this figure. The parameters 
of calculations using SBDART radiative transfer model were: 


e the water vapour column: 2.085 g/m?, 
e the ozone column: 0.35 atm-cm, 
* the tropospheric ozone: 0.0346 atm-cm, 
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Fig. 3.7 The solar flux at the snow surface (courtesy: M. Lamare) 
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Table 3.2 The coefficients of ө by Eq. (3.123) 


























PO pi À1, microns À2, microns 
3238e + 1 —1.6014033e + 5 7.95953e 4- 3 8-534e—4 ቺት- 

e the aerosol model: rural (Shettle and Fen 1979), 

9 the vertical optical depth of boundary layer at 550 nm: 0.1, 

e the altitude: 825 m, 

* the snow albedo at the surface calculated using the assumption that the grain 


diameter is equal to 0.25 mm and grains are of spherical shape. 


Taking into account that the downward solar flux appears in the nominator and 
the dominator of Eq. (3.119), the smale-scale oscillations of the curve for the calcu- 
lation of the flux is of a secondary importance. Therefore, for the back-of-e pe 
calculations one may use the approximation of the downward flux given Mio 


F(A) = pot pe "^ + pe™™, (3.123) 


where the oscillations of the downward flux are ignored and the coefficients derived 
at SZA — 60 degrees are presented in Table 3.2. 

The limits A; and A» in Eq. (3.122) can vary. In particular, CMP21 (Kipp& Zonen, 
Delft, Netherlands) pyranometers measure; 


e near IR BBA (A, = 708 nm, А = 2800 nm); 
e shortwave BBA (A, = 305 nm, А = 2400 nm). 


The visible BBA (A; = 305 nm, А = 708 nm) is obtained by subtracting the 
near IR BBA from the shortwave BBA. = 

The spherical BBA can be derived substituting 1448) into ፻9ሚ=119). In 
particular, it follows: 


ry = (А) = exp(-/a(4*)L), (3.124) 


where we used the mean value theorem for the evaluation of the integral in the 
nominator. Clearly, the wavelength A* depends on the value of L. The similar result 
holds for the polluted snow as well. 

Accurate estimations of the BBA require the numerical solution of the radiative 
transfer equation (see, e.g., Aoki et al. (2011), who also proposed the account for the 
partial illumination of a snow surface by diffuse and direct light). In particular, Aoki 
et al. (2011) has proposed the following expression for the shortwave BBA: 


. Foisfb vis ተ FNIRFb,NIR (3.125) 
b,sw — | | 
ud Fois + FwiR 
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where Ё, are global solar radiation in the visible, Fy;q is the global solar radiation 
in the NIR, 


Fb vis = (1— Justo + Juss (3.126) 
rb NIR = (1— DAS + Ріки, (3.127) 


[4 (fwin) is the diffuse fraction of the global solar radiation in the visible 
Б), is ће visible plane ВВА, ns is the visible spherical BBA, М4 9 is the 
NIR plane BBA, гу; * is the NIR spherical BBA. Sometimes the spherical (plane) 
albedo is referenced as white (black) sky albedo. 

The broadband albedo of snow is of importance for climate studies and also for 
general circulation models. It is measured at many points worldwide including the 
PROgram for Monitoring of the Greenland ICE Sheet (PROMICE) network of 
weather stations. The locations of stations are shown by red points in Fig. 3.8. The 
BBA is measured at all illumination conditions (clean, cloudy sky, see Figs. 3.8, 3.9). 


1: solar and infrared radiation, 2: tilt sensor, 3: satellite 
mms wind speed & direction, 5: snow/ice surface 
height, : air temperature & humidity, 7: ice ablation 

ose 


se’, 8; solar panel, 9: data logger, barometer and 
GPS, [0: battery, | |: ice temperature profile (8 levels) 





Fig. 3.8 The PROMICE network of weather stations in Greenland (www.promice.org) 
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Fig. 3.9 The automatic measurements of snow BBA at the KAN О PROMICE station (April 12, 
2013) (courtesy: J. P. Box) 


3.5 Finite Optically Thick Snow Layers: Reflection 
and Transmission 


3.5.1 Ambartsumian Approximation 


3.5.1.1 General Equations 


Let us assume now that the snow optical thickness (SOT) is such that snowpack 
can not be treated as a semi-infinite layer. However, it is supposed that SOT is large 
(above 5—10). Then simple approximations for light reflection and transmission by 
a snowpack can be derived. Let us show it using the Ambartsumian approximation 
(Ambartsumian 1944, 2011) for the diffuse transmittance t and the diffuse reflectance 
(spherical albedo) r of a finite turbid layer. To simplify, we study the radiative transfer 
along a line with the optical thickness x (an idealized 1-D medium). This medium 
reflects some fraction r / of initial intensity Io and transfers a fraction //6 along а 
line. 

Now we assume that this 1-D medium is composed of two 1-D media: one 
with optical thickness ту and yet another with the optical thickness t2. Clearly, the 
transmitted J, and reflected T, light intensities will given by the following equations: 


84 3 Radiative Transfer in Snowpack 


Fig. 3.10 The radiative | 
transfer in 1-1) medium —- 


I, =t(t + t2)lo, I, = r(t + т) 10. (3.128) 


We shall also introduce light intensities at the boundary between two media (see 
Fig. 3.10) and call them T, (for the radiation in the direction of incidence after 
propagating the layer of optical thickness ту) and J, (for the radiation propagating 
in the direction opposite to the incidence direction after reflection from the layer of 
optical thickness t2). 

Clearly, it follows: 


L =r(tm)h (3.129) 


and, on the other hand, 
hL =t(t)lotr(t)h. (3.130) 
Then one derives from Eqs. (3.129), (3.130): 
I; = t(1)lo +r(t)r(tm)h. (3.131) 


This means that 


/(፲1) 


| лэш 


(3.132) 


ቧበ1 


_ (to) 
І, = 1...) (3.133) 


and, therefore, 


ttt) 
/(፻] + T2) = I- rapra) (3.134) 


For the reflected radiation, we may write: 


I, = /(፲1)/2 + (т) 0. (3.135) 
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Therefore, it follows: 


г(тә)ї?(т\) 


r(t + ፻2) = (т) + 1) 


(3.136) 


The derived Eqs. (3.134), (3.136) show how the reflectance and transmittance of a 
combined layer can be expressed via reflectance and transmittance of two individual 
layers. These two sub-layers could have not only different optical thicknesses but 
also they can have different absorption and scattering properties. In a similar way, 
we may consider an arbitrary number of layers. In the discussion given below we 
assume that both media have the same scattering and absorption properties. 


3.5.1.2 Nonabsorbing Media 


Let assume that there is no absorption in the medium. Then it follows: r (т)-Е1(т) = 1 
and, therefore, 


Е t()t (12) 
“2፡5 ም 





or taking inverse values: 
ያ“ ተጄ)=ያ 0 ተያ” ሜነ – 1. (3.138) 
We may also write this equation in the following form: 
(т +m 117 0) 1+7 (о) – 1. (3.139) 
Let us introduce the function 
F(t) 21^!(r) — 1. (3.140) 
One can see that this function has the following property: 


Е(т + т) = F(11) + F(t) (3.141) 


and, therefore, 


F(t) = ст, (3.142) 
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where c does not depend on r. Clearly, this parameter must depend on the peculiarities 
of light scattering (isotropic/anisotropic light scattering, etc.) in the layer under study. 
An important point is that we arrive to the following equation: 


7! 6)= 1 ተረኛ (3.143) 


and, therefore, it follows for the transmission and reflection coefficients of a turbid 
layer: 


ct 





t(t) = e г(т) (3.144) 


+ ст = Iper 


The numerical calculations using the radiative transfer equation and also asymp- 
totic analysis of the RTE (Sobolev 1975; Minin 1988; van de Hulst 1980) show that 
the inverse transmittance can be actually presented as 


t !'(፣) = b+ cr, (3.145) 


wher{D)72 and c = 0.75(1 — g) independently on the actual phase function. 
Therefore, Eqs. (3.144) must be substituted by more accurate equations: 


CT 





1 
(тучы) 


= ; (3.146) 
+ ст b + ст 


These equations аге of great importance for the snow optics because they show 
how the diffuse reflection and transmission of a nonabsorbing snow layer (say, in the 
visible for clean snow) changes with optical thickness (assuming black underlying 
surface), if the semi-infinite layer approximation does not work. 


3.5.1.3 Absorbing Media 
Now we consider the case, when one can not neglect light absorption processes in 
an optically thick but a finite layer. Again we concentrate on the 1-D case. Let us 
assume that the incident light intensity is /o. Then the element of optical thickness 
dt absorbs the following fraction of light: 

а = (1 — оо) оат. (3.147) 


The reflected light fraction will be: 


r = (1— f)oolodr, (3.148) 
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where f describes the probability of scattering in a given direction. It is equal to f in 
the forward scattering direction and (1 — f) in the backward scattering direction for 
1-D case. Clearly, we should have due to the conservation of energy law: 


፣ተያተፀ= 1. (3.149) 
Therefore, it follows for the transmission: 
t = 1 — (1 — oof)dr. (3.150) 
Now we apply equations given above at t; = т, dt = ፻2: 


(1— f)eot2(z)dz 
1—r(t)( — odt’ 


/(፲)(1 — (1 — w f)dt 
1--፤(፻)(1 — f)ogdt ` 


r(t+dt)=r(t)+ 





(3.151) 


(т + ат) = 





(3.152) 


Let us use the linear approximation for differences dr(t) = r(t + dt) — r(t) 
and dt(t) = t(t + ат) — 1(т) with respect to dt. Then it follows: 


dr(t) _ 


(1 — fot (z), (3.153) 
dt 





dt(1) 
ч” (1 — /)шог(т)ї(т) — (1 — vof )t(1). (3.154) 





The last equation takes the following form at wọ = 1 : 





009 d - APG), (3.155) 

dt 
with a solution given by Eq. (3.144) at c = 1 — f. One can see that the derivative of 
the diffuse transmittance at c = 1 is proportional to the value of t? and not to t as 
it the case for the direct light. 

The system of two simultaneous ordinary differential Eqs. (3.153), (3.154) can be 
solved analytically. This makes it possible to derive analytical equations both for r 
and t. As a matter of fact, the system can be reduced to a single differential equation 
taking into account that the value of £? (т) in the first equation can be presented as 


2(т) = 1 + r2(z) — br, (3.156) 


where 
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_ 2(1— of) 


= I 3.157 
ሠወ9(1፦ f) dnb 


This equation can be derived by division of the first equation by the second one 
and subsequent integration. Namely, it follows: 


b 
tdt = rdr — 2ሠ (3.158) 


0፲ 


t =r? —br+C, (3.159) 
where the integration constant С = 1 because t = 1 andr = 0 at т = 0. Therefore, 
it follows from Eq. (3.153): 


dr(t) 
edt 





ያ") + br(t) — 1= 0, (3.160) 


where 
в = (1 — f)oo. (3.161) 
This nonlinear differential equation can be solved analytically. Let us introduce 


the new coordinate: z = (1 — f)wot. Then the equation to be solved can be presented 
in the following form: 





d 
0 — (r(z) _ rg )ሆርጾ) ሸዱ ro), (3.162) 


where гу! and ro are the roots of the polynomial 
VG) = 1- 2 (0) — brG). (3.163) 


This differential equation has the following solution: 


1-- - 
peg EM (3.164) 
1 — rp exp(—vz) 


where v = (rg. = ro). Let us return to the initial coordinate v. Then we derive: 


1 — exp(—2kr) 
= ፻0 2 ; 
1 — rg exp(—2kr) 





(3.165) 
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where 


1 
k= 50 = fi! — ron. (3.166) 
We derive for the transmittance from Eq. (3.116): 


1 — rë —k 
je сЕ o (3.167) 
1 — rg exp(—2kt) 





One can see from JH. 122) that ro is equal to the value of r for a semi-infinite 
layer. Therefore, it follows that rọ = e ” as discussed above. Let us introduce the 
parameter x = kt orx = y L, where L is the snow geometrical thickness and y = kk; 
is the asymptotic flux attenuation coefficient. The coefficient k is called the diffusion 
exponent. 

Then we derive: 

exp(—y)(l—exp(—2kr) ,  (l—exp(-2y) exp(—x) 


r= t= (3.168) 
1 — exp(—2kt — 2y) 1 — exp(—2y — 2x) 








or alternatively 


h 
መ с. р (3.169) 
sh(x + у) 
shy 
mM ss, 3.170 
sh(x 4- y) ( ) 


Assuming that c — 1, one derives from Eq. (3.166): 


рели (3.171) 
3(1 — g) 


The asymptotic analysis of RTE shows that in reality [see Eq. (3.62)] 
k= J3(1 — g)(1 — ox), (3.172) 


asco — 1, which means that the empirical parameter e = 1 — f given in approximate 
equations presented above must be substituted by є = 0.75(1 — g). One can see 
that reflection and transmission characteristics of a turbid layer with a given optical 
thickness for weakly absorbing media depend on just two parameters (k, y — 4 s) as 
it was in the case of asymptotical regime inside turbid layers far from its boundaries. 
The first parameter (k) describes the rate of attenuation of light field in deep layers 
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spherical albedo 





0 20 40 60 80 100 120 140 160 180 
snow geometrical thickness, mm 


Fig.3.11 The dependence of the spherical albedo on the snow geometrical thickness at the 
wavelength 450 nm [symbols—measurements (Perovich 2007), line—Eq. (3.169)] 


of turbid media. The second parameter (y) describes the light reflection from a semi- 
infinite layer of a turbid medium. Various turbid media of a given optical thickness 
t having similar values (k, y) also have similar optical properties. Equation (3.144) 
follows from Eqs. (3.169), (3.170) at probability of photon absorption equal to zero. 
More accurate Eq. (3.146) can be derived from Eqs. (3.169), (3.170), if the value of 
y in the dominator of these equations is substituted by the parameter Y — by. 

The spherical albedo r calculated using Eq. (3.169) is shown in Fig. 3.11 as 
function of the snow geometrical thickness L. Also we show the measured spherical 
albedo for the natural snow cover over various thicknesses (Perovich 2007) on the 
same figure. The value of y in Eq. (3.169) has been derived as 


y= (=). (3.173) 
Too 


where ra; is the measured spherical albedo at the largest snow geometrical thickness. 
The asymptotic flux attenuation coefficient has been derived using the following 
equation: 





m =" ге (3.174) 
— n А 
oe ia er |. 
which follows from Eq. (3.169). Here r* is the spherical albedo at the snow geomet- 
rical thickness L taken to be equal 15 mm for the experimental data shown in Fig. 3.11. 


3.5 Finite Optically Thick Snow Layers: Reflection and Transmission 91 
The value of x in Eq. (3.169) has been calculated as: 
x — yL. (3.175) 


One can see that Eq. (3.169) can indeed be used to describe the experimentally 
measured values of spherical albedo for finite snow layers. It has appeared that the 
value of AFEC is equal to 0.1464 cm! for the case shown in Fig. 3.11. The e-folding 
scale is defined as ¢ = 1/y. Therefore, it follows: ¢ = 6.8 cm. The value of ¢ gives the 
distance, where the diffuse radiation attenuates in e—times for a given snowpack. 
The transport optical thickness van be derived from the values of x = y L and y. 
Namely it follows: 


4x 
t=. (3.176) 
3y 
The transport extinction coefficient is derived as: Ку, = T/L. Also one can 


determine the transport length: /, = 1/o;,. The transport length is the length over 
which the direction of propagation of photon is randomized (see Fig. 3.12). The 
mean free path is determined by the inverse value of the extinction coefficient „уь = 
1/kex;. Then it follows: lj, = lmfp/(1 — 8). It follows that for the isotropic scattering 
both characteristic lengths (lm fps 1) coincide. For anisotropic scattering media such 
as snow /,, > уь. We shall assume that the snow asymmetry parameter in the visible 
and near infrared is equal to % as for crystalline clouds. Then it follows: /;, = 41р. 
Therefore, the direction of propagation is randomized along the distance equal to the 
4 mean free path lengths in the snowpack. 

Also the snow optical thickness can be derived assuming that the asymmetry 
parameter is known. Assuming that g — 3/4, one derives: 


И 16х 


pa ጨጨ (3.177) 
3y 
Fig. 3.12 The transport Я 
length in a highly scattering Q O 
medium 





transport length 
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The transport extinction coefficient can be derived as k,, = “=. It has appeared 
that k,- = 3.09 cm! and, therefore, extinction coefficient kext = ky, /(1— 2) is equal 
to 12.4 cnr for the case shown in Fig. 3.11. The extinction coefficient is related to 


the effective diameter of snow grains der: 





(3.178) 


where c = p;/ pj, p, is the snow density, p; is the ice density. Taking into unt that 
= 0.16 gcm^? for the experiment shown in Fig. 3.11, we derive from 3.167): 
0.42 mm, which is close to estimates of the size grains given by Perovich 
(2007). This confirms the fact that the theory presented here can be used to describe 
the spherical albedo of finite snow layers. It follows from Fig. 3.11 that the limit 
of a semi-infinite snow is reached at the depth of approximately equal to 15 cm at 
450 nm for the polluted snow studied by Perovich (2007). This means that adding 
additional snow layers under 15 cm thick snow layer will not change snow reflectance. 
Therefore, this thickness can be considered as the radiation penetration depth for the 
case studied. Although the contribution of the radiation from the depths larger than 
5 cm is already small (see Fig. 3.11). Perovich (2007) has stated that the snow with 
the depth 8 cm over the measurement platform (painted black) appears as white as 
the surrounding deeper snow. 

The theory presented above makes it possible to derive many important snow char- 
acteristics such as AFEC, transport extinction coefficient, extinction coefficient, snow 
optical thickness, and the snow grain size. The snow absorption optical thickness 
Tabs = Text L can be also derived: 

xy 
Tabs = 4 ° (3.179) 

We show the spherical albedo calculated at various values of diffusion exponent k 
in Fig. 3.13. The experimental results are also given there. The snow optical thickness 
can be presented as 


= (3.180) 
t= E w 
a ЕЛ 
where 
dg (3.181) 
Di 


pi — 0.917 вст? is the bulk ice density and p; is the snow density. The value of p; 
was equal to 0.16 ест! for the experiment shown in Fig. 3.11. Therefore, it follows: 
A — 0.522. We have also estimated the ice grain diameter from the measurements as 
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spherical albedo 





0 20 40 60 80 100 120 140 160 180 200 
optical thickness 


Fig. 3.13 The dependence of the spherical albedo on the snow optical thickness at the wavelength 
450 nm and various values of the diffusion constant k [symbols—measurements (Perovich 2007), 
lines—Eq. (3.169)]. The fitting (see the dashed line) experimental results at two points shown by 
circles leads to the estimated value of the diffusion exponent equal to 0.012 


discussed above. It appears that the effective grain diameter is eaual to 0.42 mm. 
Therefore, it follows: т = çL, where the scaling constant ç = 1.23 Ç Dre have used this 
scaling constant to present experimental measurements of the spherical albedo as the 
function of snow optical thickness in Fig. 3.13. The deviation of experimental data 
from the dashed line shown in Fig. 3.11 is probably due to vertical variation of snow 
properties and also due to the errors of the measurements (albedo, snow geometrical 
thickness). It should be pointed that the snow optical thickness at which the limit of 
a semi-infinite snow is reached depends on the diffusion exponent k being smaller 
at larger values of k. One can state that the semi-infinite medium case is reached at 
SOT — 120 for the experiment shown in Fig. 3.11. This limiting value of SOT is 
much smaller in the near IR part of the electromagnetic spectrum. 


3.5.2 Reflection and Transmission Functions 
of Nonabsorbing Snow Layers 


Let us consider now the transmission function for an optically thick snow layer. 
Because ፲ > 1, the dependence of transmittance on the azimuth can be neglected. 
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Also both theory and observations show that the angular distribution of the trans- 
mitted light can be presented as a linear function of the cosine of the observation 
angle. Itis possible to show that the normalized angular distribution of the transmitted 
light (to its value in zenith direction) does not depend on the value of t. Taking into 
account that we consider homogeneous layers, the principle of reciprocity must be 
satisfied. This means that 


To(5, n, т) = to(r)uo($ )uo(n). (3.182) 


The separation of the angular and geometrical coordinates is one of main results 
valid for optically thick scattering layers. 
The diffuse transmittance is defined as (for any value of single scattering albedo) 


t= ‚С (3.183) 
0 


L= 
tP (£) = 24 TS: n)ndn, (3.184) 
where 
= ] 2x 
T(E, n) = 2x / TE, n, ф)аф. (3.185) 


Taking into account the considerations presented in the previous section, we 
derive: 


1 
То, n, t) = b+ pa, በን (3.186) 


where b = 1.072, р = 0.75 and т, = kl is so-called snow scaled or transport 
optical thickness, / is the geometrical thickness of snow layer and kp = (1 — g)kex; 
is the transport extinction coefficient. This equation can be used to study transmission 
function of a finite nonabsorbing snow layer. To derive k;,, one needs to measure the 
transmittance and snow geometrical thickness L: 

Ttr 


Or = T> (3.187) 


where 


4 
Tir = zm ® n, T)uo(E)Juo(n) 53 (3.188) 
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and we assume black underlying surface. Similar equations can be derived if not the 
transmission function but rather the plane transmittance f, (£) or diffuse transmittance 
t are measured. 

The accuracy of Eq. (3.186) as compared to the solution of RTE is given in 
Fig. 3.14. It follows that the transmission function depends on a single parameter- 
scaled optical thickness till ዌራ = 1. The linear dependence of 1/T on the scaled 
optical thickness is valid till т,. = 1.5 for the case studied. 

Let us consider now the reflection function R(E, 7, т) assuming that there is no 
lightabsorption in the medium. Approximate equation for this function can be derived 
in the following way. 

Let us take the semi-infinite nonabsorbing medium and make an artificial cut 
inside the medium at large optical thickness t. Then the reflection function of a 
semi-infinite layer can be presented as the sum of the reflection function of optically 
thick upper layer and also the transmission function of radiation diffusely reflected 
from the lower semi-infinite nonabsorbing layer. Taking onto account that the diffuse 
reflection coefficient of the lower layer is 1, we derive: 


Reoo (Š, n, p, т) = RoE, n, т) + To, n, т). (3.189) 


- -- - approximation, g=0.85 
exact result, g=0.85 
exact result, g=0.75 

zenith observations 

SZA=60 degrees 


inverse transmission function 





0 1 2 3 4 5 6 7 8 9 10 
scaled optical thickness 


Fig. 3.14 The dependence of the inverse value of the transmission function on the scaled optical 
thickness calculated using Eqs. (3.186) and the numerical RTE solution at single scattering albedo 
wo = 1, the Henyey-Greenstein phase function and two values of the asymmetry parameter (0.75, 
0.85) at the zenith direction. The solar zenith angle is equal to 60°. 
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This means that 


፳0(፳, n, т) = Roo (Š, n, 9) — To(E, n, т). (3.190) 


Therefore, the calculation of R(é, n, т) for the case of nonabsorbing medium is 
reduced to the calculation of the functions studied above. Because the dependence of 
Roxoo(&, n, @, т) on the phase function can be neglected, one can use the parameter- 
ization to determine R..o(&, 7, ф) for a given geometry and then derive 7o(&, n, т) 
and transport optical thickness from the reflectance measurements of a given finite 
snow layer. 


3.5.3 Reflection and Transmission Functions of Weakly 
Absorbing Snow Layers 


In the case of weakly absorbing layers, one can use similar equation as given ለ ው 
However, the diffuse transmittance coefficient must be substituted by that presentet 
Then it follows (Zege et al. 1991): 


T(&, n, т) ио(&)ио(2), (3.191) 


Е shy 

- sh(x + by) 
where we have accounted for the fact that b differs from unity. The equation for 
the reflectance of a finite weakly absorbing layer can be derived in a similar way as 


discussed above. Then it follows: 


R. (5, n, 9) = RE, п, т) + TE, п, z) exp(—x — y). (3.192) 


Therefore, one derives: 





exp(—x — у)ио(&)ио(1), 
(3.193) 


shy 
RE, n, Фф, z) = ፳ኋ0(. n, 9) exp(—uy) PENT 


where u = ио(&)ио(т) Rao (£, n, ф) and we have accounted for «አዲ... We see 
that the problem of finding R(E, n, o, т) isreducedto the calculation ofthe previously 
studied functions R.o(&, n, 9) and uo (ë). 
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3.5.4 The Optically Thick Snow Layers with Arbitrary Level 
of Absorption 


The equations for the arbitrary level of absorption can be presented in the following 
form (Sobolev 1975, 1984): 





T = а 3.194 
E, n, т) = Г Neip(-2kp) GA. (3.194) 
R(E, n, 9, 1) = Ro ፦ T(E, n)N exp(—x), (3.195) 
where 

1 
M= 2 | Ponnan, (3.196) 

--1 

1 

N= 2 f исуї=т)лат. (3.197) 


0 


The escape function u (m) at any level of light absorption in the medium has been 
studied in detail by Yanovitskij (1997). It can be derived if the intensities in deep 
layers i (5) and also the azimuthally averaged reflection functions of a semi-infinite 
turbid layer are known. 

It follows that the simple relationships given above (at b — 1) can be derived from 
more general expressions given in this section under assumptions: 


M = 1 — №, N = exp( y). u(&)u(n) = uo(Ë)uo(n), (3.198) 


which are valid as oo — 1. 


3.5.5 Account for Underlying Surface and Vertical Snow 
Inhomogeneity 


Equations for finite layers given above can be applied for black underlying surface. 
An account for Lambertian surface with albedo A can be done in the following way 
(Liou 1992, 2002): 


ATE) 


R = R ack , * 
(£, ቫ› т) black (Š, Nn, t) ተ ied 


(3.199) 
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Here r is the spherical albedo of a snow layer for the case of black underlying 
surface, (ድ) is total (direct + diffuse transmittance of a snow layer for the case 
of black underlying surface), Ry. (5, 7, т) is the snow reflection function for the 
case of black underlying surface and we have assumed that a snow layer is vertically 
homogeneous. Therefore, R(E, 7, т) can be presented via kn elationships valid 
for the case of black underlying surface. The derivation of E 148) is simple. In 
particular, multiple reflections from the Lambertian surface can be presented in the 
following form (Stokes 1862): 


RE, n, т) = Rusa (Ë, n, т) + £6)AQ + Ar + (Ағ)? +... Jt). (3.200) 


The summation of this series gives the result presented in Eq. ሟቅ In a similar 
way one derives for the transmission function: 


Ar Et) 


T ‚1, = Тыас EE 
(E, n, т) = Tua (Ë, n, т) + icm 


(3.201) 


Snow is usually composed of sub-layers with various concentration of impurities 
and snow properties such as density, ice grain shapes and sizes. In the absence of 
absorption, the vertical inhomogeneity can be treated using the average values of the 
snow phase function/asymmetry parameter. However, if absorption must be taken into 
account, one must find the solution of the radiative transfer equation with vertically 
varying phase function, single scattering albedo, and extinction coefficient. The task 
must be solved using the numerical calculations. 

Alternatively, the vertically inhomogeneous medium can be substituted by a 
homogeneous one with a special choice of effective single scattering albedo (Yanovit- 
skii 1991). The problem of two-layered snow with upper optically thick layer and 
lower semi-infinite layer with albedo гоо can be considered using equations given 
above. In particular, it follows in this case: 


К(&, n, v) = R, (£, n, t) + Toota (E )tu (n) (3.202) 


1 — тогы 


where index u signifies the upper layer. The account for the light scattering and 
absorption effects in atmosphere overlying the snowpack under assumption of 
the Lambertian reflective ground can be done with Eqs. (3.199), (3.201) as well. 
In this specific case the first terms in Eqs. (3.199), (3.201) must be substituted by 
the atmospheric reflection and transmission functions, respectively, and A has the 
meaning of the snow surface albedo. The plane albedo гу, total transmittance t, and 
spherical albedo r must be changed to these parameters for the overlying atmosphere 
(Kokhanovsky et al. 2005b). 

The correct treatment of the vertically inhomogeneous snow layers can be 
achieved by the numerical solution of the RTE for given vertical distribution of 
local optical characteristics of snow layers such as single scattering albedo, phase 
function, and extinction coefficient. 
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3.6 The Polarization of Light Reflected from Snow 


The light reflected from the snow surface is polarized. It is characterized not only by 
its intensity / but also by other characteristics such as the degree of polarization P, the 
ellipticity e, and the orientation of the polarization plane x. These characteristics can 


be derived from the Stokes-vector—parameter T with the components / , Q,U,V: 


УО? +U2+ v? 1 U 1 V 
P= сий шл ‚у= 3 wcan( ©) e=- "|. |) (3.203) 


ያ 2 ያ 





The Stokes—vector components can be expressed via the components of the 


electric vector É perpendicular (E,) and parallel (Ej) to a selected plane (Stokes 
1852): 


I = ЕЕ + E*E,, Q = ЕЕ — E EU = ETE, + ЕЕ, 
V =i(E*E, — E*E,). (3.204) 


The ellipticity of reflected light is small and can be neglected for clean snow 
layers. Although the ellipticity of reflected light can be different from zero due to 
the presence of living cells in snow. 

For most of applications, one can assume that light reflected from the snow 
is partially linearly polarized. The degree of polarization is small in the visible. 
However, it increases in the ice absorption bands. The example of the measured 
degree of linear polarization of light reflected from a semi-infinite snow layer is 
given in Fig. 3.15 for the solar zenith angle in the range 55—58 degrees. The illu- 
mination from the Sun comes from the backward direction (see Fig. 3.15) along 
the principal plane (PP), which contains the perpendicular to the surface and the 
direction towards the Sun (Tanikawa et al. 2014). 

The degree of linear polarization is defined as 


P = —, 3.205 
7 ( ) 


where Y = y Q? + U? is the polarized reflectance. One can also introduce the 
bidirectional polarized reflectance factor (BPRF) П = л Y/& Fo, where л Fo is the 
incident solar flux at the unit area perpendicular to the light beam and š is the cosine 
of the solar zenith angle. 

The parameter U is close to zero in the principal plane (PP) containing the surface 
normal and the Sun direction (Peltoniemi et al. 2009). Therefore, the degree of linear 
polarization shown in Fig. 3.15 is derived from the measurements as follows: 


P2--, 2 
: (3.206) 
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Fig. 3.15 Viewing angle dependence of the degree of linear polarization at five selected wave- 
lengths. Positive viewing angles correspond to the forward direction and the negative ones corre- 
spond to the backward direction. The arrow on each figure shows the position of the Sun. The 
illumination from the Sun comes from the backward direction along the principal plane, which 
contains the perpendicular to the surface and the direction towards the Sun (Tanikawa et al. 2014) 


The positive values of P correspond to the light polarization in the plane perpen- 
dicular to the principal plane (у = 7/2). The negative values correspond to the 
polarization in the plane parallel to the scattering plane (y = 0). It follows from 
experimental data presented in Fig. 3.15 that light reflected from snow is weakly 
polarized in the visible. The degree of polarization increases in the near infrared, 
where light is polarized in the plane perpendicular to the principal plane for most of 
geometries. Although there are angular regions, where light is polarized in the plane 
parallel to the principal plane. It follows from Fig. 3.15 that there are neutral points, 
where light reflected form snow is completely unpolarized (P — 0). The lower degree 
of polarization in the visible as compared to the near infrared can be explained by 
the fact that the absorption of ice is small in the visible and multiple scattering domi- 
nates. Also light polarization by irregularly shaped ice crystals is weak (see Fig. 2-10). 
Multiple scattering leads to the increase of entropy and decrease of degree of polar- 
ization. Therefore, we conclude that spectral regions, where snow appear brighter 
are characterized by lower levels of degree of polarization. The peak at the viewing 
zenith angle around 55 degree (glint region) is due to the fact that the solar zenith 
angle is approximately equal to the Brewster angle o = arctan(n), where n is ice 
refractive index, for the experiment presented. The melt-freeze crust is characterized 
by the largest values of the degree of polarization. It is a thin, glittering ice layer, 
which is sometimes visible on the snow on sunny days. Therefore, the polarimetric 
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measurements at the Brewster geometry can be used to classify the snow conditions 
using remote sensing ground and satellite instrumentation. 

The experimental measurements of light polarization by snow layers can be under- 
stood in terms of vector radiative transfer theory. In particular, it is known that the 
Stokes vector of light reflected and transmitted by snow layers is satisfied to the 
following vector radiative transfer equation (VRTE) (Liou 1992): 


>> 

a (Š Р , "ev 

00 = -& T Gb» | 5ጻ.(8.8)67ወ jd. (3.207) 
An 


The matrix £., transforms the Stokes vector defined in the meridional plane holding 


the normal to the scattering layer and the direction d to the Stokes vector in the 
scattering plane. The matrix Lo is needed for the transformation of the Stokes vector 
of scattered light from the scattering plane back to the meridional plane (see Appendix 
2). This is due to the fact that the matrix ĉsca in Eq. (3-194) is defined in the scattering 
plane and the vector I is defined with respect to the meridional plane. For randomly 
oriented ice grains the matrix 6,,, is reduced to a scalar value k,,, and 


dT (8) 


ሙያ T@+ /፡ሯጮ [s] Эё (Q ya Š. (3.208) 


4л 


We will be mostly concerned with solutions of Eq. (3.195) for a plane-parallel 
snow layer illuminated by the Sun. The interaction of solar radiation with extended 
snow fields is well-approximated by the solution of this idealized problem. The solar 
light with the zenith angle 90 uniformly illuminates a plane-parallel scattering layer 
from above. We will assume that properties ofthe layer do not change in the horizontal 
direction. Then properties of light field depend only on the vertical coordinate Z 
and the direction б, specified by the zenith angle 9 and the azimuth o. The main 
task is to calculate angular distributions d (9, o, z). Usually only measurements of 
Y (9, o, 0) at the top of the snow (reflected light) are performed. Therefore, we will 
be concerned with the angular distribution of the reflected light. 

Equation (3.195) takes the following form for the plane-parallel snow layer: 


d I (9. 
cos 2 KELON = —Ooxt T (9, ф) 


dz 
Ju В. v, ф 9, Ф) ÈI, ф'). (3.209) 


The most simple case to study is that of idealized homogeneous snow layers with 
values of kext and бу not dependent on the position inside the ctoud. Then we have 
from Eq. (3.198): 
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d Í (ә, . . 
eg oo -T (8,9) + efu 126. 9, ф) ЁГ", Ф), 








(3.210) 
where т = kextz is the optical depth, 
^ AMO seal’, o > 9, 
P= | ? 9) (3.211) 
Куса 
is the phase matrix and 
Куса (3 212) 
оо = ; 
ን 


is the single scattering albedo. The phase matrix, extinction coefficient and single 
scattering albedo of snow grains have been studied in Chap. 2. 

It is useful to distinguish between diffused is and direct (or coherent) Т. == 
Fo (и — uo) (o — qo) light in the general solution Т (0, ф). Неге Fa describes 
the Stokes vector of the incident attenuated flux. It is assumed that the layer is 
illuminated in the direction defined by the incidence zenith angle #6 = arccos(uo) 
and the azimuthal angle #0. The multiply scattered light is observed in the direction 
specified by the zenith observation angle 9 = агссоѕ(и) and the azimuth o. Namely, 
we write: 


— — — 
I (0,9) = 140,4) + 1.0, 9). (3.213) 


The substitution of this formula in Eq. (3.199)-gives 


OL P(o, g > 8,9)L Tawo’) 





T 
+7 — B, (Oo, Qo > v, ф) 7,መ[-ሯ=ኗ) (3.214) 
cos Yo 


Here P, is the correspondent scattering matrix defined with respect to the merid- 


ional plane, Т (Fo, Qo, Uo, Vo) is the Stokes vector of incident light flux at the top 
of a cloud. It follows for unpolarized solar light: Qo — Uo — Vo. The solution of this 
equation under boundary conditions stating that there is no diffused light entering the 


cloud from above and below, allows to find Т a (9, ф). Also it follows that (ә, ф) 
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is given simply by: 





=> => 
I .(0, ф) = F (0, @)ó(cos 9 — costo) (ф — qo) ew(- ) (3.215) 


COS 0 


The solution of Eq. (3.213) is a more easy task than that of Eq. (3.199) because 
we avoid the necessity to deal with the divergence in the direction of incident light 
existing in the framework of Eq. (3.199). There are multiple numerical algorithms, 
which can be used to solve the integro—differential VRTE (3.213) accounting also 
for the vertical variation of the snow properties (http://en.wikipedia.org/wiki/Atmosp 
heric radiative transfer codes). The input for the corresponding codes in the case of 
vertically homogeneous snow layers is the optical thickness, single scattering albedo 
and the phase matrix of snow grains, which depend on the size of grains, their shapes 
and wavelength (Siewert 2000). 

Equations for optically thick layers presented in the previous section can be gener- 
alized to account for polarization. Correspondent equations for azimuthally averaged 
reflection R and transmission 7 matrices were obtained by Domke (1978a, b). They 
have the following forms for isotropic symmetric light scattering media: 


RE, n) = R< (£, n) — NT (E, n) exp(—kz), (3.216) 


M exp(—kr) 


Fajr 3217 
ушы | 6) (3.217) 


T(E,n) = 





where only two-dimensional matrices and vectors are involved. Other components 
of generally four-dimensional matrices and vectors vanish due the azimuthal aver- 
aging. Note that this is also the case for a normal illumination of a snow layer. Then 
the azimuth does not enter theory at all. Here Roo (€, n) is the azimuthally averaged 
reflection matrix of a semi-infinite medium with the same optical characteristics as 
a finite slab under study. One can also introduce the vector 6 (u), which describes 
the intensity and degree of light polarization in deep layers of a semi-infinite scat- 
tering medium [in a so-called asymptotic regime, when the intensity and polarization 
angular distributions are symmetrical with respect to the normal to a scattering layer 
and exponentially decrease with the depth (~ exp(—K7))]. Functions Roo (£, n) and 
Ө (2) determine all parameters in the equations given above. The explicit equations 
for the calculation of Roo (£, ዛ) and O(n) and also the escape vector Е (п) аге 
discussed by Domke (1978а, b, 2016). 

We see, therefore, that intensity and polarization characteristics of reflected and 
transmitted light for optically thick snow layers are determined by the reflection 
matrix of a semi-infinite layer and the angular distribution of the light intensity and 
polarization in deep layers of the same medium. This reduction of a problem for a 
finite optical thick snow layer to the case of a semi-infinite medium is of a general 


^ > 
importance for the radiative transfer theory. Note that the matrix Rə and vector © 
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are obtained from solutions of the well-known integral equations, which can be found 
elsewhere (van de Hulst 1920—4e Rooij 1985). 

Equations (3.215), (3.21 valid only for the azimuthally averaged matrices. In 
practice, however, measurements are performed for a fixed azimuth. The transmission 
matrix is azimuthally independent in the case of optically thick layers. The azimuthal 
dependence in the reflected light disappears in some specific cases (e.g., for the case 
of normal illumination of a now surface). 

Equations 3.215),(3.21 С} simplified for nonabsorbing media. Then it follows: 





RE, n) = RE, 8) — TE, n), (3.218) 
A 4 — —T 
TE, n) = IESER Ko) Ko (n), (3.219) 
where 
1 
2 22T > 

qo = pem duy Ko (n) j. (3.220) 

0 

Here 


; [1 
J= (3) (3.221) 


is the unity vector, g is the asymmetry parameter, 


1 
— 


3 А = 
Кой) =F) n2 | dE’ (ё, n) |j (3.222) 
0 


and RO. (u, шо) is the azimuthally averaged reflection matrix of a semi-infinite nonab- 
sorbing medium. This matrix is completely determined by the phase matrix Ê, intro- 
duced above. It does not depend on the single scattering albedo and optical thickness 
by definition. Clearly, the first component of the vector Ko coincides with the escape 
function uo discussed above. 

These asymptotic equations are simple in form. However, they can be used only 
if auxiliary functions and parameters are known. Their calculations, however, can 
be quite a complex procedure. However, it appears that for weakly absorbing media, 
when single scattering albedo is close to one, simplifications are possible. Then it 
follows (Kokhanovsky 2003): 


RGE, п) = RooolE, M exp(-yÂE, m) – FE, mMer- у, (3223) 
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76. 8) =t RoE) E, m, (3.224) 


where x = kt, y = 319, К = 430-w- а), DE, n) = 


5-1 ጭ =»ታ E sinh y : : В 
Ry. n) KoE) Ko Q), = ааваа у 15 the global transmittance of a scattering 


1 
layer, g = 5 f uo] n^ dy ~ 1.07, and Ro (t, Ho) is the reflection matrix of a semi- 
0 


infinite nonabsorbing layer with the same phase matrix as an absorbing layer of a 


finite thickness under study. The two-dimensional vector Kolu) describes the polar- 
ization and intensity of light in the Milne problem for nonabsorbing semi-infinite 
media (Wauben 1992). The components Ко (ш) апа Ко (и) of this vector were 
calculated by Chandrasekhar (1950) for Rayeigh particles (g — 0) and by Wauben 
(1992) for spherical les with the refractive index п = 1.44 and the gamma 
particle size distributi: 5) wit 11.3, ag = 0.83 jum. The wavelength of À 
was equal to 0.55 jum. Note, that t del of spheres with и = 11.3, ғо = 0.83 um, 
and n — 1.44 is generally used to characterize particles in clouds on Venus (Hansen 
and Travis 1974). It follows for the effective size а, у, the effective variance A,;, and 
the asymmetry parameter g, respectively, in this case: Asf = 1.05ит, Aef = 0.07, 
g — 0.718. It was found that the ratio p; — — 20 which gives the degree of polar- 
ization for transmitted light is very low. It changes from zero to 1.2% while the 
escape angle changes from 0 till 90 degrees. Note that for Rayleigh scattering we 
have a change from 0 till 11.7% for the same conditions. This means that light trans- 
mitted by thick snow layers is almost unpolarized. It is possible to understand this 
on general grounds. Indeed, the polarization of unpolarized solar light occurs due to 
single scattering events. Multiple light scattering leads to an increase of entropy and 
the reduction of initial polarization arising in single scattering events. 
Equations given above can be simplified for nonabsorbing media (y — 0): 


RE, n) = Ros (Е, m) — T(E, n), (3.225) 
^ => —T 
T(é,n) =t Ko(€) Ko Q), (3.226) 
where 
t= РЕ. Е (3.227) 
mm" 21(1 — g) i 


is the global transmittance. 

Let us apply Eq. (3.225) to a particular problem, namely, to the derivation of a 
relation between the spherical albedo r = 1 — t and the degree of polarization of 
reflected light p; (u) at the illumination along the normal to the snow layer (uo = 1) 
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by a wide, unidirectional unpolarized light beam. The value of p; (m) is given simply 
by — R2 (1, ш)/ Ru, ш) in this case. Thus, it follows from Eq. (3.225): 


pico (Š ) 
= —— — __, 3.228 
”/(፳) I-ü- 086) ( ) 
where 

uo(1)uos (u) 
N ED 3.229 
OY" ж are 

Roo A, 
Pi») = a (3.230) 


and we accounted for the equality: Ko2(1) = 0. It appears that the value of N(ë) is 
close to | for most of observation angles, which implies the inverse proportionality 
between the brightness of a turbid medium and the degree of polarization of reflected 
light (rp; ^ pico). This inverse proportionality between the spherical albedo r and 
the degree of polarization p; was discovered experimentally by Umow (1905). Equa- 
tion (3.228) can be considered as a manifestation of this important law, which has 
important applications in reflectance spectroscopy (Hapke 1993). 

Equation (3.228) is easily generalized to account for the absorption of light in a 
medium using the exponential approximation described above. Namely, it follows: 


Pino (Š) 





р(&) = [= ከቦር የመ (3.231) 

where 
TS M 3.232 
^ sinh(x + ey) ሂን 

and 
።0(1)#0(፳) 

N*(é) = ———*. 3.233 
(ë) RECE) ( ) 


Values of př (£) and ፳፲ (1, £) represent the degree of polarization and reflection 
function of a semi-infinite weakly absorbing medium at the nadir illumination. Note, 
that Eq. (3.317) can be written in the following form: 


PIE) = c(Š, T) Pio (ë), (3.234) 


where 
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1 
c@ v) = ТЕ ao (3.235) 





can be interpreted as the polarization enhancement factor, which is solely due to 
a finite snow depth. It follows for semi-infinite layers that the transmittance f is 
equal to zero and c = 1 as it should be. Also it follows from Eq. (3.320) that 
zeroes of polarization curves for semi-infinite and optically thick finite snow layers 
almost coincide. This is due to the fact that the function N*(£) only weakly depends 
on the angle. Multiple light scattering fails to produce the polarization of incident 
unpolarized light. It only diminishes the polarization of singly scattered light. Thus, 
the angles where polarization is equal to zero for semi-infinite layers are almost equal 
to those for the case of single light scattering. 

In the case of semi-infinite weakly absorbing snow layers, it follows 
(Kokhanovsky 2003): 


RoE, m) = Res (E, т) exp(—yDE, m) (3.236) 
or assuming that y — 0: 


А ^ — —T 
RE, n) = Ros (Š, n) — y Ко(&) K ( (n). (3.237) 


Let us write Eq. (3.237) in the explicit form: 


en | _ Gee | _ ‚о 10 
፳ኋዕ21 Roo22 Rooc21 R0oc22 Ko2(&)Koi(n) Ko (£) Koo (n) J 
(3.238) 


In the case of unpolarized solar light illumination we need to multiple Eq. (3.238) 
by the column vector (van de Hulst 1980) 


+ fà 
= Er (3.239) 
Then it follows: 
Ru Коп Ко (&) Кот (п) 
= = 3.240 
E» bol i) 
or 


Root = Коо — y Koi (E) Koi (n), (3.241) 
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Коо = Rooo21 — Y Koo (6) Koi (m). (3.242) 


Equation (2.241) gives the relationship between the reflection function of 
absorbing (Ro; = R5511) and nonabsorbing (Rog = Rooo11) snow layers as у 0. 
The second equation shows the difference of reflection functions for the polarizations 
parallel (R;) and perpendicular (R,) to the principal plane: 


Roo2t = Root — Roo, = Коо — y Koo (E) Koi m). (3.243) 


Also it follows for the degree of polarization: 


Poo = Pooo + Ay, (3.244) 
where 
K K K K 
a = KOKO p | ROKO (3248) 
Roco Roco 


Taking into account that К (£) describes the degree of polarization in deep layers 
of nonabsorbing media, we can assume that K2(Ẹ) << 1 and derive from Eqs. 
(2.244), (2.245): 


Poo — m. (1 + жк) (3.246) 
Ооо 
or 
Poo = Poco exp) J (3.247) 


It follows from Eq. (2.247) that the degree of polarization increases with the level 
of light absorption in the medium. This is confirmed by the experimental data shown 
in Fig. 3.15. It should be pointed out that the approximation of spherical particles 
can not be used in snow polarimetry. The peak in the degree of polarization due to 
enhanced light polarization in the rainbow region common for large spherical parti- 
cles does not exist for snow surfaces. The dependence of the degree of polarization 
on the viewing zenith angle calculated using Eq. (2-247) is presented in Fig. 3.16. 
The results presented in Fig. 3.16 closely reproduce experimentally derived features 
of DOLP given in Fig. 3.15 such as low values of degree of polarization for snow, 
the decrease of degree of polarization towards the antisolar point (-57 degrees) and 
neutral point around -40 degrees for the SZA = 57 degrees. 

The experimentally measured and calculated using Eq. (2-247) spectral degree of 
polarization is shown in Fig. 3.17. The value of pos; has been calculated using single 
scattering approximation. For the nadir illumination case as shown in Fig. 3.17 the 
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Fig. 3.16 The degree of polarization calculated using Eq. (2.247) at the SZA = 57 degrees and the 
fractal model of ice grains (Macke et al. 1996) under assumption that y = (eL, where о is the 
bulk ice absorption coefficient and L — 1 mm. The wavelengths are 0.55 (lower curve), 1.02, 1.24, 
and 1.6 jum (upper curve) with larger wavelengths corresponding to the larger values of degree of 
polarization. The calculations are performed in the principal plane with negative viewing angles 
corresponding to the backscattering direction 


value of poss coincides with the normalized phase matrix element Pj?(see Appendix 
2) for nonabsorbing ice crystals in the approximation under study. We have used the 
geometrical optics calculations of Р performed by Macke et al. (1996) for fractal 
ice grains. 

It follows that Eq. (2.247) provides means for the calculation of the snow degree 
of polarization at channels, where the snow can be considered as a weakly absorbing 
medium (spherical albedo above 0.1). The theoretical results around 2000 nm are 
not reliable due to the fact that snow can be not be considered as a weakly absorbing 
medium at this channel. Taking into account that it follows for the reflection function: 


Ki) KG) | (3.248) 


Roo = Ross apl =y R 
Ооо 


we conclude that the second Stokes parameter 


Q = — poo Ræ = — poco Roco (3.249) 
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Fig. 3.17 The spectral degree of polarization calculated using Eq. (2.247) at the УХА = 50 and 60 
degrees at the nadir illumination in the principal plane (forward scattering) and the fractal model of 
ice grains (Macke et 81. 1996) under assumption that у = VoL, where a is the bulk ice absorption 
coefficient and L = 1mm. The results of experimental measurements (Sun et al. 2021) are also 
shown 


does not depend on the single scattering albedo for the wavelengths, where snow is 
weakly absorbing. This is also confirmed by experimental data (Sun et al. 2021). 


It follows from Eq. (3.249): 


Poo = ዥ”› (3.250) 


where А = рого Roo. Equation (3.250) is the manifestation of ће Umow law, which 
state that the degree of polarization is inversely proportional to the brightness of the 
turbid medium. 
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Chapter 4 (R) 
Remote Sensing of Snow EM 


4.1 Determination of Local Optical Parameters of Snow 


4.1.1 Semi-Infinite Snow Layers 


Spectroscopy is used to refer to the measurement of transmitted/reflected radiation 
intensity as a function of wavelength. The intensity of light transmitted/reflected 
by snow layers depends on local optical parameters of snow such as extinction 
coefficient, absorption coefficient, and directional light scattering coefficient (phase 
function). In turn these local optical parameters depend on the snow microstructure 
and pollution load. Therefore, several important snow parameters can be derived 
from snow spectral reflectance and transmittance measurements. In some cases the 
spectral polarization parameters are measured. Usually, the snow reflectance spec- 
troscopy is used due to much easier experimental setup as compared to the snow 
transmittance/internal light filed measurements. However, diffuse light transmit- 
tance measurements enhance the information content and, therefore, provide addi- 
tional information for the solution of inverse problems of snow optics. Therefore, 
we consider the techniques for the determination of snow properties based on both 
reflectance and transmittance measurements. 

We shall start from the determination of snow local optical characteristics using 
spectral reflectance and transmittance measurements. Snow microstructure param- 
eters and pollution load can be derived, if snow spectral optical characteristics 
are known. The main optical characteristics of snow are the extinction coefficient, 
absorption coefficient, and phase function. 

Snow is composed of ice crystals. Therefore, the snow phase function is close to 
that of crystalline clouds. Phase functions of various crystalline clouds have been 
studied experimentally. They can be modelled by the combination of exponentials 
or alternatively by the combination of the Henyey-Greenstein phase functions as 
discussed in Chap. 2. Clearly, the morphology of snow covers is not identical to that 
of crystalline clouds due the snow metamorphism processes and external factors such 
as wind, ground temperature variations, solar insolation, and precipitation including 
© Springer Nature Switzerland AG 2021 115 
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rain. However, one may expect that phase functions of snow are close to that if 
ice clouds with irregularly shaped crystals measured in the laboratory. The main 
feature of snow phase function is its smooth behavior in the backward region. The 
asymmetry parameter of ice clouds has been measured experimentally and it appears 
that its value is close to 0.75 in the visible. This means that the average geometro- 
optical cosine of scattering angle is close to 0.5 for nonabsorbing snow. The value 
of g increases for absorbing grains. It is close to 0.97 for strongly absorbing convex 
ice particles in random orientation. 

The extinction coefficient of snow is difficult to measure directly. However, again 
one may expect that extinction coefficient of light in snow can be modelled in the 
same way as for ice clouds. Then it follows as discussed above: 


3Cice 
def 





; (4.1) 


kext = 


where d, is the effective diameter of grains and Cice is the volumetric concentration 
of snow crystals. Taking into account that often Cice ^: 1/3, we find that and snow 
optical thickness (SOT) т ~ L/d,y, where L is the geometrical thickness of the snow 
layer. Therefore, one may expect that just 5—10 layers of snow particles constitute 
optically thick medium with the optical thickness above 5-10. This is confirmed also 
by the naked eye observations of thin snow layers on, say, black underlying surfaces. 
Usually just several cm of snow is enough to make underlying surface invisible. 

It follows from the discussion presented in the previous Chapter that radiative 
transfer characteristics of weakly absorbing snow layers are determined by mainly 
two parameters, namely: similarity parameter s and the asymptotic flux attenuation 
coefficient y. The similarity parameter s can be derived from the measurements of 
the diffuse snow reflectance r under diffuse illumination conditions in the visible 
and near-infrared as discussed above: 


1 
#= “ገ (4.2) 


A similar equation can be derived, if the plane albedo is measured. The value of 
y is obtained from the measurements of the light flux in snow (say, at two levels 
in snow, see Eq. 3.74). The spectra $ (А), y (À) determine the spectral behaviour of 
snow spectral reflectance and transmittance in the visible and near-infrared. 

It should be pointed out that several local optical characteristics of snow are 
related at small values of probability of photon absorption. For instance, it follows: 


ፅ — k В Каьѕ Kabs y 
s= = =” = = = = (4.3) 
30—8) 30-8) k КК y 3k,, 


y = kkext, k = y3B(1— 8), (4.4) 
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kir = (1 e g)kexi. (4.5) 


This means that the snow transport extinction coefficient can be derived from the 
pair (s, y): 





ACP) 
К (A) = А 4.6 
r) = z 6) (4.6) 
Also it follows for the snow absorption coefficient: 
ka О) = у О)5 (À). (4.7) 


4.1.2 Finite Snow Layers 


4.1.2.1 Reflectance and Transmittance Measurements 


Ithas been demonstrated above how the snow optical parameters can be derived from 
light reflection measurements (for a semi-infinite snow layer) in combination with 
studies of light fluxes in deep snow layers. A similar information can be derived from 
the measurements of simultaneous spectral transmittance and reflection of optically 
thick finite snow layers of a given snow depth. Let us assume that the diffuse reflection 
andtransmission coefficients of snow are measured. They are defined by the following 
equations for the snow layer with the geometrical thickness L as discussed in the 
previous Chapter: 


shx shy 


— , t, = ——..U (4.8) 
sh(x + by) sh(x + by) 


r |, = 


where x = y L, y = 4s. To simplify the analytical solution of the inverse problem 
we ignore the difference of the parameter. m unity. The pair (x, y) can be derived 
from (r, t) given by Eq. (4.8) at la ically. Let us show it. First of all, we 
note that 


t h 
LII (4.9) 
r shx 

where we removed the subscripts to simplify the equations. Also it follows that 


sh(x + y) 2 shxchy + chxshy (4.10) 


and, therefore, 


(shx chy + chx shy)t = shy, (4.11) 
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(спу shy + VI+ shx shy)r = shy, (4.12) 
rchy +ty1+sh?x = 1, (4.13) 

rchy + М? — r2 + r2ch2y = 1, (4.12) 
chy = = (4.15) 


Taking into account that 


hx = “shy, (4.16) 


we derive: 





] + r2 — 232 
Ges t (4.17) 
t 4r? 





Therefore, it follows finally: 


1 2 — ሥን) 1-22 2 
x= Е — ) Е = ы; (4.18) 








Alternatively, we can write (Хере et al. 1980): 


=] dde. 161... Еа (419) 
х=; ap |? "| 3p 42 | : 





where 

b= /(1 + r2 — ይ) — 4r2 (4.20) 
and we used the equality: 

arshA = ክ(4 + у 1 + A2). (4.21) 


We can derive the pair (y, s): 


(4.22) 


and 
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4 
kabs = Yy/4, ky = £. (4.23) 
3y 


4.1.2.2 Reflectance Measurements of Snow At Two Values 
of Geometrical Thickness 


The pair (y, y) can be also derived from just reflectance or transmittance measure- 
ments at two (or more) values of the geometrical thickness of a snow layer. Let us 
show it. 

The snow spherical albedo can be presented in the following form: 


Fs = lg — t exp(—x — y) (4.24) 


where 


‚— (17 APC) exp(=x)_ (4.25) 
1 — exp(—2y — 2x) 





Then measuring the snow reflectance for a semi-infinite snow layer (defined as a 
layer which does not change with the further increase of snow thickness), we derive 


y = In(r;?). (4.26) 


The dependence of the parameter у = 4 s on r; is given in Fig. 4.1. 
Performing measurements of the snow reflectance (with black underlying surface 
at a fixed snow geometrical thickness h), we derive: 


1፦ r (rs 
2yL)= Z = 38 
exp( yL) 1 — r, (A)Prso0 





(4.27) 


and, therefore, 





| In 77! 


= 4.28 
y 2L (4.28) 


If it is needed, the measurements can be performed at two arbitrary geometrical 
thicknesses (e.g., not for a semi-infinite snow). Then the pair (y, y) can be also 
derived. 
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Fig. 4.1 The dependence of the retrieved value of y on the spherical albedo r, 


4.1.2.3 Transmittance Measurements of Snow At Two Values 
of Geometrical Thickness 


The parameters (y, y) can be also derived from snow transmittance measurements 
at two snow thicknesses A1, h2: 


mu (1 — exp(—2y)) ехр(—х1) P (1 — exp(—2y)) ехр(—хо) (4.29) 
1 — exp(—2y — 2x1) 1 — exp(—2y — 232) 








Namely, it follows: 


1—пе% — Da he 
e2 = ш i (4.30) 
1-በሪ”፤ 1--ርሮ= 





It means that 


(1— te") — he’) — (1— 5e?)(1—tje 1) = 0 (4.31) 


or, therefore, 


ty osh(y (E — Iz) — tish(yz) + 5А (ух) = 0, (4.32) 
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where z = #1, £ = h>/ hi. This transcendent equation can be used to derive у. Then 
it follows: 
ыы ае (4.33) 
= — In| ——_ |, s 
) 3 eges 


Assuming that £ = 2, one derives: 





th — tj + 2t5ch(yz) = 0 (4.34) 
or 
= E (4.35) 
z 
where 
t — tt 
=> » x (4.36) 


4.1.2.4 Nonabsorbing Snow Layers 


In the visible, clean snow is nonabsorbing and, therefore, kaps = 0 and k; can be 
derived from the transmittance measurements at the snow sample of thickness h. 
Namely, it follows: 


EI 
kir = ar Uo (ና, n.B)uo(E)uo(n) — а), (4.37) 


where Tọ ME n, , B) is the inverse snow transmission function. 
Let us assume that the snow diffuse transmission coefficient t (£, h) is measured. 
Then it follows: 


ky = ae, h)ug(§) — e). (4.38) 


The equation is even simpler if the diffuse transmission under diffuse illumination 
(t) is measured: 


ky = 7 (h 4.39 
акті (й) — е). (4.39) 
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4.0 Snow Grain Size Retrieval 


Snow is composed of irregularly shaped ice crystals of various forms. As a matter of 
fact, each crystal is unique and has its own specific shape and geometrical dimensions. 
Therefore, the notion of snow grain size is not well defined for snow. For the size of 
particles, one may choose, e.g., the average maximal/minimal dimension of crystals, 
which can be measured using a lens. Also more advanced approaches to characterize 
crystals are available. Reflectance spectroscopy makes it possible to determine the 
effective absorption length (EAL) in snow. This parameter can be used to characterize 
the snow by optical methods. In particular, it follows for the spherical albedo of a 
semi-infinite snow layers 85 wọ — 1: 


diss exp — eQ), (4.40) 


where о (À) is the bulk ice absorption coefficient and £ is the effective absorption 
length. Опе can see that £ can be easily derived from the measurements of r, (А) at 
a particular wavelength A: 


_ In?r, 
- ዐዐ) 





(441) 


This equation is valid only for clean snow and must be used in the near-infrared, 
where r, deviates from 1 (say, at 1020 nm) and the theory used to derive this equation 
remains valid. The value of £ is an important parameter characterizing any snowpack. 
It can be related to the size of grains using the following equation: 


dey = б (4.42) 


_ 9-8) 
S = 168 


where (see Eqs. (3.119), (3.120)) 


(4.43) 


The value of £ can be also found from plane albedo measurements: 


t= omn _ (4.44) 
(Aus (Ho) 


and also from the reflectance measurements: 


_ In?(R/Ro) 


= Gar” (4.45) 
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where Ro is the nonabsorbing snow reflectance (in the visible). For a vertically homo- 
geneous snow, the derived parameters £ and d, shall not depend on the wavelength. 
However, due to different penetration of radiation with different wavelength to snow- 
pack, both 4, у and £ depend on the wavelength in case of vertically inhomogeneous 
layers. This feature can be used to study the vertical profiles of correspondent param- 
eters. For better understanding of snow vertical inhomogeneity and buried polluted 
snow layers, one needs to study the snow profiles along the vertical snow walls. 

The snow grain diameter can be also retrieved using measurements not at a single 
wavelength but also in the ice absorption bands centered at specific wavelengths. 
The spectral observations of the snow surface can be used to assess the vertical 
distribution of snow grains in the snow top layer owing to the different penetration 
length of radiation depending on the wavelength. 


4.3 Determination of Snow Specific Surface Area 


Specific surface area is the property of solids defined as the total surface area of 
material per unit of mass. Therefore, it follows in the case of snowpack: 
NS 
y = E (4.46) 
NM 





where N is the number of ice grains in the unit volume, S is their average 
surface area and 


M = рУ (4.47) 


is their average mass, p; = 0.917 g/cm? is ће ice density and V is the average volume 
of grains. Therefore, it follows for the snow specific surface area: 


S 
gu =, (4.48) 
piV 
Let us introduce the effective grain diameter using the following equation: 
ү 
def = E (4.49) 


In the case of monodispersed spheres this diameter coincides with the actual 
diameter of spherical particles. Then it follows: 





(4.50) 
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Therefore, one can see that the snow specific surface area can be derived if the 
average grain diameter is known and vice versa. As it has been discussed above, 
the measurement of the effective grain size for snow is a difficult matter. Therefore, 
usually direct measurement techniques are used to estimate SSA. In particular, the 
methane adsorption (proportional to 5) technique is used to estimate S. The mass of 
a snow sample can be found by weighting. Then one derives SSA as the ratio shown 
in Eq. (4.46). 

This is quite an involved procedure, therefore, it looks as an attractive solution of 
the problem to develop the technique for the determination of SSA using the snow 
effective absorption length derived from snow reflectance measurements. It has been 
demonstrated that £ can be presented in the following form: 


2BV 
= Е (4.51) 
3(1 — g)® 


where Ф is the average surface area of grains, g is the asymmetry parameter, and B 
is the absorption enhancement parameter. The accuracy of this expression is high for 
weakly absorbing snow (for the wavelengths below 1000 nm or so). It also follows: 


32B 
L= —— (4.52) 
3(1 — g)pio 
Or 
32B 
o= —__.. (4.53) 
3(1 — g)pil 
Equation (4.53) can be presented in the following general form: 
A 
c= —, (4.54) 
pil 
where 
32B 
А = —— (4.55) 
3d = 8) 
depends on the snow type. It follows from Eq. (4.54): 
Ino = € — lnh, (4.56) 


where 
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h = pit (4.57) 


and 


e=InA. (4.58) 


It is essential to study the value of A for different types of snow experimentally. 


4.4 Determination of Snow Impurity Content 


4.4.1 General Equations 


Snow originates from the atmosphere, where various aerosol (solid and water-soluble) 
particles are present. Therefore, it contains not only frozen water but also other 
chemical substances such as brown/black carbon, dust, litter from the trees, etc. Also 
various types of algae in various quantities depending on location can be present. 
Reflectance spectroscopy is traditionally used to understand the chemical composi- 
tion of various substances. In this section we consider the application of reflectance 
spectroscopy to derive the chemical composition and concentration of various impu- 
rities in snow. The variation of pollution load in various snow samples can be assessed 
from Table 4.1. 

Let us assume that the diffuse reflectance of snow under diffuse illumination 
conditions for semi-infinite is studied. Such a scheme makes it possible to avoid the 
use of auxiliary functions, which are potential sources of an error. The measured 
spherical albedo can be presented as 


rs = exp(—y). (4.59) 


Therefore, the spectral values of the similarity parameter y can be directly retrieved 
from the measured reflectance spectrum. In particular, it follows: 


yQ) = In(r; 'ዕ.))› (4.60) 


Therefore, the inverse problem involving multiple light scattering in snow is 
reduced to the determination ofthe type of pollutants and their load from the similarity 


parameter 
yA) = 4 | TERT (4.61) 
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Table 4.1 „Тһе шелше d Impurity Concentration Location Source 
concentration of impurities at 
several locations Soot 0.1—0.3 ng/g South Pole Flanner et al. 
(2007) 
Soot 1-30 ng/g Summit, Flanner et al. 
Greenland (2007) 
Soot 28-210 ng/g Rural Flanner et al. 
Michigan (2007) 
Soot 80-826 ng/g French Alps | Flanner et al. 
(2007) 
Soot 17—5700 ng/g Urban Flanner et al. 
Michigan (2007) 
Dust 0.04—0.1 mg/g European Biagio et al. 
Alps (2015) 
Dust 20—50 mg/g Alpine and Skiles et al. 
subalpine (2012) 
Algae 0.1-1.5 cells/mL | Antarctic Gray et al. 
penisula (2020) 
Algae 1-25cells/mL Greenland Stibal et al. 
ice sheet (2015) 











More precisely, the spectral probability of photon absorption 6 must be studied 
because for most of cases (but not always!) light scattering in snow is governed 
mostly by ice grains and not by impurities. Therefore, we can derive from Eq. (4.61) 


ВО) = O1n’(r,(A)), (4.62) 


where the parameter Ө = 3(1 — g)/4 can be derived assuming the asymmetry 
parameter for ice grains. The error in the estimated © leads to the uncertainties in 
the derived PPA and also concentration of impurities. In particular, we may assume 
that g = 3/4 in the visible and near-infrared and, therefore, Ө = 3 

Now we need to relate the PPA A (À) to the concentration of po ts. This can 
be done in the following way. Let us assume that pollutants are externally mixed 
with snow grains and their influences on scattering (and light extinction) in snow is 
negligible as compared to the contribution of snow grains. Then it follows: 


kii A) + kans 0) 
pe сн, 


ext 


(4.63) 


where 


i 3Cice ice im 
Геи = ШЕ" Obs (A) = Bcijceotice (A), К (A) = Cpol K pol (A). (4.64) 
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Here, Cice is the volumetric concentration of ice in snow, Сри is the volumetric 
concentration of pollutants in snow, Фу is the effective grain diameter, о (A) is the 
bulk ice absorption coefficient, B is the absorption enhancement factor, and K po; (A) 
is the volumetric absorption coefficient of pollutants: 

— pol 
K poi (A) = “ሠ. (4.65) 
pol 


spol . : А ud 
where p is the average absorption cross section of pollutants and V м is the 


average volume of impurity particles. One may also introduce mass absorption 
coefficient of pollutants: 


K ро 
K pol,mass (A) = P E. (4.66) 
D pol 





where оро is the density of a pollutant. Therefore, it follows for the probability of 
photon absorption: 


det 
BQ) = [Bai Q) + CK poi (A) ፻፪ (4.67) 


where ር = Cpoi/Cice is the relative concentration of pollutants. One can see that the 
value of c can be derived from PPA: 





1 
c= E — Ba. | (4.68) 
det K poi 
or 
3B 
c= ar p 4.69 
K poide f ( ) 


if the wavelength selected for the inversion procedure is in the visible, where one 
can neglect the second term in Eq. (4.68) due to small light absorption by ice in the 
visible. One can see that the relative concentration of pollutants can be derived if the 
parameters В, d, f , and K por (or alternatively, K pot,mass) are known. The pair (В, dey) 
in Eq. (4.69) can be derived as discussed above. The main problem in the precise 
determination of c is the determination of the volumetric absorption coefficient of 
pollutants (see Eq. 4.69). The value of K po; depends on the size distribution f (a) of 
pollutant particles and also on their complex refractive index m pot = при — i X pol- 
This makes retrievals of c problematic because both f(a) and m po: (and, therefore, 
K pot) are not known in advance. Let us consider this problem for main snow pollutants 
such as soot and dust. 
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4.4.2 Soot 


In the case of soot, the scatterers are much smaller than the wavelength of incident 
light and, therefore, 





እ aga (A) = WO pol (А)У pots (4.70) 
4л X pol 
K poi (A) = Чоро (A), A pot (A) = UW (4.71) 
where 
9 DO. 
y = id . (4.72) 


2 
2 2 2 2 
("ы +1— 35 + 467 oi em 


There is some uncertainty with respect to the spectral refractive index of soot. 
Assuming that пры = 1.75, Хра = 0.45, we derive that Y = 0.9. Finally, it follows 
from Eq. (4.69) for the relative concentration of soot in snow: 


ፅ 


c — II , 
а pou deg 





(4.73) 


where П = 3/ У. Asa matter of fact, some of soot can be internally mixed. Therefore, 
the concentration as derived from this equation can be biased. 


4.4.3 Dust 


The derivation of the pollution load for the case of dust is even more complicated 

because: (1) dust can influence both light scattering and absorption processes in 

snow; (2) the spectral refractive index of dust differs depending on the origin of dust. 
The spherical albedo of dust-loaded snow can be presented as 





ANY 
r(a) = exp] —é — + cpako( =) IL , (4.74) 
0 


where l, = 1/(1 — g)kex; is the transport extinction length, € = 4/ МЗ, ko is the 
volumetric absorption coefficient of impurities at the wavelength Ао. 
It follows from Eq. (4.74): 


4.4 Determination of Snow Impurity Content 129 


r(À) = exp} — (so +9[፲) ) ; (4.75) 





where 
2 ko 
í = CiceBE lj, а = p? (4.76) 


ር = Cpol Сісе 18 the relative concentration of pollutants. It follows that the spectrum 
of the dust loaded snow depends on three parameters: v, q, £. They can be derived 
from the measurements of spherical albedo at three wavelengths. Namely it follows 
from Eq. (4.75): 








v= Gy" == a). (4.77) 
ln A Aicel A3 
where 
012 AL " _ Inr) 
р = Іп со (2) а= ағ)" (4.78) 


Equation (4.77) can be used to derive ће snow grain size of polluted snow and 
also the concentration of pollutants c. Namely, it follows from Eq. (4.76): 


9(1— 2) B 
ie PT 4.79 
16В "c Fc pte) 





One can assume that ር = 9 and В = 1.6 (Kokhanovsky et al. 20218). This 
makes it possible to estimate the snow grain size from the measurements of EAL. 
For the determination of the concentration of pollutants we need to know the value 
of ko. Clearly, the value of ko is correlated with the value of the absorption Angstrom 
parameter v. The following correlation equation derived from the Mie theory under 
certain assumptions with respect to the size distribution of dust grains and also their 
complex refractive index can be used (Kokhanovsky et al. 20218): 


2 
ko = >` b; v, (4.80) 
i=0 


where bo = 10.916 -2.0831, b2 = 0.5441 and the volumetric absorption 
coefficient is express 1/mm. 

The accuracy of the theory developed above is demonstrated in Fig. 4.1, where the 
following triplet of the wavelength has been used to retrieve v, q, £ from the measured 
plane albedo spectra: 410, 500, and 825 nm. One can see that three-parameter formula 
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(4.75) can be used to determine the snow plane albedo 
pep (4.81) 


The various parameters of snow derived from the spectra shown in Fig. 4.2 are 
presented in Table 4.2. The relative mass concentration is found as follows: 


ንም. (482) 


i 


where py = 2.65g/cm? is the density of dust and p; = 0.917g/cm? is the density 
of ice. 
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Fig. 4.2 The spectral snow albedo: measurements-red lines, calculations using Eq. (24)—blue 
lines. Dates refer to the acquisitions of the albedo spectra (Kokhanovsky et al. 20218) 


Table 4.2 The parameters of dust-loaded snow derived at three dates at the Torgnon site 
(Kokhanovsky et al. 2021a) 





SZA, degrees 
16.05.2018 | 24.44 3.00 |2.391e-5 
17.05.2018 | 27.21 2.51 | 1.5176-4 
18.05.2018. 26.98 3.36 |2.304e-4 


q, l/mm 





9.63 1.15 11.7 
9.11 1.60 71.4 
10.11 2.33 106.9 
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The formulae given above are valid only for external mixture of pollutants and 
ice grains. In practice, part of particles can be externally mixed and another part of 
impurity particles can be contained inside (internal mixture) ice grains (increasing 
absorption), which complicates the retrieval procedure. In addition, several and not 
just one pollutant (including algae) may present in snow making retrievals even more 
involved. 


4.5 Spaceborne Remote Sensing of Snow 


4.5.1 Spaceborne Instrumentation 


Snow covers huge polar regions and also it has considerable extent in the Northern 
hemisphere during winter. Snow permanently presents almost over whole Antarc- 
tica, which is the fifth-largest continent almost twice of Australia (14,200,000 km?). 
The same is true for Greenland (2,166,086 km7). Therefore, it is of importance 
to develop techniques for snow monitoring using satellite instrumentation. Usually 
passive optical (visible and near — infrared) and microwave sensors are used for snow 
monitoring from space. In particular, snow water equivalent can be derived from 
passive microwave observations. In this chapter we consider only the optical snow 
remote sensing. The list of optical instruments most often used for snow monitoring 
is given in Table 4.3. They are also briefly discussed below. 

Multiple Spectral Imager (MSI oard Sentinel-2 operates in a broad spectral 
range 0.443-2.2 рт providing me ments at a high spatial resolution (10—60 m, 
depending on the channel). The 12 detectors on each focal plane are mounted in a 
staggered formation to cover the whole 20.6? instrument field of view, resulting in a 
compound swath width of 290 km on the ground track. 

Ocean and Land Color Instrument (OLCI) is a single-view push-broom imaging 
spectrometer, which measures the top-of-atmosphere reflected light in the spectral 
range 400-1020 nm with a spatial resolution of 300 m. It operates on board Sentinel-3. 
The swath is 1270 km. 

Sea and Land Surface Temperature Radiometer (SLSTR) on board Sentinel-3 uses 
two independent scan mirrors rotating in opposite directions each scanning at a rate 
of 200 scans per minute providing dual-view capability. The swath width is 1400 km 
for the nadir observation and 700 km for the along track view. The spatial resolution 
is 500 m in solar channels. It performs measurements in visible, short/mid wave 
infrared. It also has thermal channels located at 10.85 and 12.02 um. 

Second Generation Global Imager (S-GLI) is a push-broom instrument providing 
11 channels in the visible and near infrared spectral region and 6 channels in the 
shortwave infrared and thermal infrared regions of electromagnetic spectrum. Two 
channels have a 3-view capability and perform measurements of the polarization 
state of reflected light, which is crucial for atmospheric correction procedures. The 
spatial resolution is 0.25-1 km depending on the channel. 
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Visible Infrared Imaging Radiometer Suite (VIIRS) observes and collects global 
satellite observations that span the visible and infrared wavelengths across land, 
ocean, and atmosphere. A whiskbroom radiometer by design, it has 22 channels 
ranging from 0.41 jum to 12.01 jum. Five of these channels are high-resolution 
image bands or L-bands, and sixteen serve as moderate-resolution bands or M- 
bands. VIIRS is one of five instruments onboard the Suomi National Polar-orbiting 
Partnership (NPP) satellite platform that was launched on October 28, 2011. SNPP 
(formerly called the National Polar-orbiting Operational Environmental Satellite 
System (NPOESS) Preparatory Project) serves as a bridge between the Earth 
Observing System (EOS) satellites and the next-generation NASA-NOAA Joint 
Polar Satellite System (JPSS). Four JPSS missions have been planned to last through 
2031, and each of them will host a VIIRS instrument as part of their payload. JPSS-1 
(NOAA-20) was launched on November 18, 2017. The swath is 3000 km with nearly 
global coverage every day (1:30p.m. local solar time). 

The Moderate Resolution Imaging Spectroradiometer (MODIS) was launched 
into Earth orbit by NASA in 1999 on board of Terra (EOS AM) and in 2002 on board 
the Aqua (EOS PM) satellites. The instruments capture data in 36 spectral bands 
from 0.4 to 14.4 microns at varying spatial resolutions (2 bands at 250 m, 5 bands 
at 500 m, and 29 bands at 1 km). Together these instruments image the entire Earth 
every 1—2 days. 








4.5.2 Cloud Screening 


The retrieval of snow properties using spaceborne observations requires the applica- 
tion of cloud screening and atmospheric correction algorithms. The cloud screening 
over dark ocean surfaces could be easily performed because clouds are much brighter 
as compared to underlying ocean surface. The snow and cloud surfaces are both 
bright. Therefore, the cloud screening over snow is a complex matter. Several tests 
are used to select the cloudless satellite pixels. They are summarized in Table 4.4. The 
tests are based on the fact that the reflected solar radiation at the top-of-atmosphere 
for scenes with snow and clouds have different spectral, directional, and polarization 
characteristics. These differences are explored to distinguish underlying snow and 
cloud surfaces. The most important feature used for cloud identification over bright 
surfaces is screening of gaseous atmosphere underneath the clouds for cloudy scenes. 
This leads to the increase of top-of-atmosphere reflectance in the absorption bands of 
gases with maxima of concentrations close to the terrestrial surface such as oxygen 
and water vapor. The threshold values (THV) used to distinguish the cloud scenes 
depend on the underlying surface height, geometry of observation, and channel. 
They can be found preparing the histogram of TOA reflectances at specific gaseous 
absorption wavelengths for clear and cloudy scenes. Such histograms show clear 
separation of two regimes (cloudy, clean), which can be used te the determination of 
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Table 4.4 Cloud tests based on the value of reflectances R at several spectral regions and observa- 
tion geometries. The value of A is the standard deviation of reflectances measured at a given point 
at different times or at a given time for an area surrounding the location of interest. The tests 6 and 


7 refer to the reflectances at the scattering angle 0 defined after Eq.(3.32) 




















N |Name Test Physical property 

1 Oxygen A-band test R(761 nm) > THV1 Larger reflectance for cloudy 
scenes in the gaseous 
absorption bands due to 
screening of gases under the 
clouds 

2 Маѓег vapor absorption test | R(940 nm) > THV2 See above 

R(1340 nm) > THV3 
3 Spectral reflectance check | R(1020 nm) < THV4 Smaller absorption by particles 
R(1640 nm) < THV5 in clouds as compared to 

crystals in snow due to their 
smaller sizes 

4 Spatial inhomogeneity A > THV6 High spatial variability of 

check cloudy scenes 
5 Temporal inhomogeneity | A > THV7 High temporal variability of 
check ground scenes covered by 

clouds 

6 | Angular reflectance check | R(148°)/R(120°) > ТНУ8 | Higher variability of angular 
distribution of reflected 
radiation for clouds (e.g., 
rainbow, glories) 

7 Polarization check R(148?)/R(120?) > THV9 | See above except for the 
degree of polarization P 














THVs and the development of corresponding cloud screening procedures. Another 
possibility is to use the TOA reflectance in the absorption bands of water and ice. 
Water droplets and ice crystals in clouds are smaller as compared to the size of snow 
grains, therefore, the respective TOA reflectances are larger for the cloudy scenes. 
The THVs can be determined in the same way as described above. Also the liquid 
water clouds can be identified by the shift of the absorption band as compared to 
the case of ice grains in snow (see Appendix). In addition water clouds have unique 
directional and polarization characteristics (rainbow, glory) not existing for snow 
fields. Terrestrial surfaces have unique features for a given location. The existence 
of these features in the satellite images points to the absence of clouds. In addition, 
one can check satellite images at a given location with frequent re-visit time. The 
cloudy pixels are characterized by larger temporal variability as compared to clean 
scenes. All features underlined above can be used to prepare a cloud mask and select 
cloud free scenes for the subsequent snow property retrievals. 
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4.5.3 Atmospheric Correction 


The optical signals detected by the satellite instrumentation are influenced not only by 
the underlying snow surface. The characteristics of the propagation channel (atmo- 
sphere) is of importance as well. Therefore, the spectral and angular distribution of 
Stokes parameters of reflected solar light at the top-of-atmosphere must be measured. 
This makes it possible to determine both atmosphere and underlying surface proper- 
ties for a given satellite scene. In most of cases the information content of measure- 
ments is low and one should use additional databases or assumptions to make the 
determination of snow characteristics from spaceborne measurements. 

Often the spectral measurements of the TOA reflectance R(A) are performed. 
Then assuming that the underlying surface is Lambertian, one can use the following 
relationship between R(A) and albedo of underlying surface ғ, (А) (Liou 1992): 


= rs(A)Ta(A) 
R=R,+ n Or 0)" (4.83) 


where 7; is the atmospheric transmittance on the way from the TOA to the surface and 
back to the satellite position, А. is the atmospheric reflectance for the case of black 
underlying surface and r, is the spherical albedo of atmosphere. These functions 
depend on the aerosol/molecular optical thickness, phase function, and atmospheric 
single scattering albedo. The spectral spherical albedo of underlying surface can be 
derived from Eq. (4.83) analytically: 


R(A) — Ra) 
Ta(A) + (RQ) — R (А) а(А) 





rà) = (4.84) 


Let us assume that the atmospheric optical thickness is equal to zero. Then it 
follows: Ra(à) = ra(A) = 0,፲;(እ) = 1, and r, = R(A) as it should be in the 
assumption of underlying Lambertian surface. If the atmospheric optical thickness 
is not zero, the functions R, (А), T; (À), andr, (A) can be estimated using the radiative 
transfer modelling (see Appendix). 


4.5.4 Snow Albedo, Snow Grain Size and Snow Specific Area 


Yet another approach for atmospheric correction is the use of near-infrared channels. 
In this case the influence of the propagation channel is often can be neglected and 
one can assume that the satellite measures directly the snow reflectance function, 
which can be presented as (see Eqs. 2-101, 3.109, 4.42): 


R = Ro EA ЫШ (4.85) 
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In particular, if one uses the channels A, = 865 and Az = 1020 nm, one can estimate 
the parameters Ro and í in general ex ion (4.85) from reflectance measurements. 
Namely, Eq. (4.85) can be used for the e of £ (at 1020 nm) and Ko = Ri Rm 46 m 
1/(1 — b), b = х/о Газ, where indices signify the wavelengths used. This makes it 
possible to determine the snow reflection function at all channels including visible 
and UV wavelengths using Eq. (4.85) and the derived parameters Ro and £. In this 
way one can determine not only spectral and broadband snow albedo but also one can 
derive snow grain size and snow specific surface area using the approach described 
for the case of ground observations at the beginning of this Chapter. The typical OLCI 
spectrum over snow is given in Fig. 4.3. The channel located at 1020 nm is sensitive 
to the size of crystals in snow. The retrieval results based on OLCI measurements 
are shown in Fig. 4.4. 
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Fig. 4.3 OLCI TOA spectrum over snow at solar zenith angle (SZA) equal to 63.61°, viewing 
zenith angle (VZA) equal to 20.63°, solar azimuthal angle (SAA) 54.14 degrees, and viewing 
azimuthal angle (VAA) of —64.25°. The relative azimuthal angle is equal to 118.39 degrees. The 
OLCI measurements have been performed over Dome C (Antarctica) on November 10, 2017. The 
regions affected by various atmospheric gases are shown (Kokhanovsky et al. 2021b) 
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Fig. 4.4 Sentinel-3 daily mosaic of a snow grain size and b broadband albedo over Greenland for 
July 28, 2017 (Kokhanovsky et al. 2019) 


In addition, the difference between measured TOA reflectance (see Fig. 4.3) and 
calculated snow reflectances in the visible and UV can be used to determine the 
aerosol optical thickness and also total ozone. Such an approach is valid only for 
clean snow surfaces. Otherwise, the spectrum of reflected light governed not only by 
the atmospheric effects and properties of snow grains. The absorption and scattering 
of light by impurities in snow must be accounted as well. 

The inter-comparison of snow broadband albedo retrieved from space and 
measured on ground in the framework SICE algorithm (Kokhanovsky et al. 2020) is 
shown in Fig. 4.5. Onecan see that the satellite observations provide similar values of 
BBA as compared to ground observations. The SICE BBA albedo retrieval procedure 
relies on the measurements of TOA reflectance at 865 and 1020 nm with subsequent 
determination of spectral and snow albedo as described above. Also the snow grain 
size and specific surface area can be retrieved from optical spaceborne measure- 
ments. The results of satellite retrievals of SSA are compared to those derived from 
ground observations are shown in Fig. 4.6. One can see close agreement of satellite- 
derived and ground-measured of SSA. The monthly grain diameter, snow specific 
surface area, albedo at 1020 nm and broadband albedo for July 2017 over Greenland 
as derived from ОСІ js shown in Fig. 4.7. 
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SCO ህ site: 72.3942 М, 27.2593 W. 988 m 
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Fig. 4.5 The broadband albedo time series from the SCO О PROMICE automatic weather station 
site in comparison to the broad band retrievals from the SICE retrieval, the OLCI retrieval after 
(Kokhanovsky et al. 2020) (SSSNOW) and for the NASA MODIS MOD10A1 product (Hall and 
Riggs 2016) 
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Fig. 4.6 The temporal trends of satellite derived (OLCI) and in situ SSA for 2017 and 2018 years 
in Greenland (75.625 N, 35.973 W) (Kokhanovsky et al. 2020) 
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Fig. 4.7 Monthly grain diameter a snow specific surface area b albedo at 1020 „A broadband 
albedo e for July 2017 as derived from OLCI (Kokhanovsky et al. 2019) 
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4.5.5 Snow Fraction and Snow Extent 


Snow extent and albedo are main characteristics of snow cover as far as climatic 
effects of snow are of concern. Snow extent is monitored on a global scale using 
satellite observations in a broad range of electromagnetic spectrum from UV to 
microwaves. Microwave instruments can be used for monitoring snow cover under 
cloud fields. Optical methods can be used for studies of snow cover for clear sky 
only. The spatial resolution of passive microwave—derived snow cover and snow 
water equivalent is rather coarse (25 km). This is not the case for optical measure- 
ments, where the observations on the scale of several hundred meters are common 
(depending on the sensor). 

Normalized differential snow index is used as a main parameter to distinguish 
snow free and snow-covered surfaces. It is defined as 

Ri — Кә 


NDSI = —— ——, (4.86) 
Ri + К 


where R; is the top-of-atmosphere (TOA) reflectance in the visible and R> is the 
TOA reflectance in the near-infrared. Depending on the remote sensing instrument 
specifications, various bands are selected. In particular, MODIS NDSI relies on 
the channels located at 0.55 and 1.6 jum. The reflectance К» is small over snow. 
Therefore, pixels with relatively large values of NDSI are assumed to be snow covered 
(in addition to the large the reflectance values of reflectance at 865 nm). Various 
threshold values are employed in various algorithms. MODIS snow extent algorithm 
relies on several spectral tests including NDSI » 0.4, R(650 nm) » 0.1, R(865 nm) » 
0.11 to detect snow in a non-densely forested region (Hall et al. 2015). 

In some cases, the ground scene is only partially covered by snow. Then one can 
assume that the reflectance R can be presented as a weighted sum of reflectances of 
snow (R,) and underlying free of snow surface (Re): 


RA) = (1 — f) RQ) + f Rs), (4.87) 


where f is the snow fraction and it is assumed that the atmospheric correction of 
the TOA reflectance R has been performed. One can see that the measured and 
atmospherically corrected reflectance is a linear function of the snow fraction f. It 
follows from Eq. (4.87): 


. RQ) - RA) 
ዚወጸ = 


Therefore, the subpixel snow fraction can be estimated, if the snow and underlying 
free of snow surface reflectances are known. The example of spectral albedos (angular 
integrals of reflectance) of various snow and underlying free of snow surfaces are 
shown in Fig. 4.8. 
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Fig. 4.8 The spectral albedo 
of soil, vegetation and snow Clean snow: 
surfaces with various 
effective grain radii re. The . г = 250 um 

case of dust-polluted snow is r,=1000 um 

also shown. The vertical ux en rt 
lines correspond to the laser Soil 
wavelengths used for the 

snow depth estimation 
(Deems et al. 2013) 


Vegetation 


Spectral albedo 





Wavelength (um) 


The variability of the clean snow reflectance R, in the visible is low. Therefore, it 
can be estimated using radiative transfer calculations. Alternatively, R, can be taken 
from the brightest pixels surrounding the pixel under study. The value of А, can be 
taken from the same pixel except recorded before snowfall. If the measurements are 
performed at the near infrared wavelength, where the snow reflectance is low (say, 
at 2.1 um), one derives from Eq. (4.87): 


Ry 
=1-።-፦. 4.89 
/ Renir ( ) 


In this case one needs to know the underlying surface reflectance before the 
snowfall Re nir- As a matter of fact, the reflectances of many surfaces in the NIR and 
visible can be often related. For instance, Kaufman and al. (2002) used the linear 
relationship 


Re nir = y Revis» (4.90) 


where y © 2. Assuming that measurements are performed in the visible (А1) and 
near infrared (A5), one derives from Eq. (4.87): 


ROA) = (1 — РК (А) + f R, (À), (4.91) 


КО») = (1 — f) RQ2) + f R, (А). (4.92) 


Assuming that R,(A2) < < 1 and using Eq. (4.90), one derives from Eqs. (4.91, 
4.92): 
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R(3) = (1 — f)y Re(A1). (4.93) 

It follows from Eqs. (4.91, 4.93): 
RQ4) = RQ3)/y + ЈЕ, (А1) (4.94) 


or 


_ ROO) — RQa)/y 


= : 4.95 
5 ፳,(፡.1) eS 


where А, (А1) can be taken from the surrounding brightest pixel, К(А1) апа R(A5) are 
measured reflectances at two wavelengths and y depends on the underlying surface 


type. 
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Appendix 


A.1 Complex Refractive Index of Ice 


In this Appendix we present the spectral complex refractive index of ice m — n-ix. 
It must be pointed out that the imaginary part of ice refractive index is very small 
in the range 300—600 nm (in the optical range). Therefore, for many purposes, one 
may assume that it is equal to zero. However, the precise value of x is needed for 
several applications including the determination of snow and ice impurity load using 
the measurements in the visible region of the electromagnetic spectrum, visible light 
























































Table А.1 The spectral. እ, pm h X 

dependence of complex ice 

refractive index m — n-ix 3.000E-001 1.3339 2.0E-011 

(Warren and Brandt 2008) 3.500E-001 1.3249 2.0E-011 
3.900E-001 1.3203 2.0E-011 
4.000E-001 1.3194 2.365E-011 
4.100E-001 1.3185 2.669E-011 
4.200E-001 1.3177 3.135E-011 
4.300E-001 1.3170 4.140E-011 
4.400E-001 1.3163 6.268E-011 
4.500E-001 1.3157 9.239E-011 
4.600E-001 1.3151 1.325E-010 
4.700E-001 1.3145 1.956E-010 
4.800E-001 1.3140 2.861E-010 
4.900E-001 1.3135 4.172E-010 
5.000E-001 1.3130 5.889E-010 
5.100E-001 1.3126 8.036E-010 
5.200E-001 1.3121 1.076E-009 

(continued) 
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Table A.1 (continued) 











































































































Appendix 

X, итп n X 

5.300E-001 1.3117 1.409E-009 
5.400E-001 1.3114 1.813E-009 
5.500E-001 1.3110 2.289E-009 
5.600E-001 1.3106 2.839E-009 
5.700E-001 1.3103 3.461E-009 
5.800E-001 1.3100 4.159E-009 
5.900E-001 1.3097 4.930E-009 
6.000E-001 1.3094 5.730E-009 
6.100E-001 1.3091 6.890E-009 
6.200E-001 1.3088 8.580E-009 
6.300E-001 1.3085 1.040E-008 
6.400E-001 1.3083 1.220E-008 
6.500E-001 1.3080 1.430E-008 
6.600E-001 1.3078 1.660E-008 
6.700E-001 1.3076 1.890E-008 
6.800E-001 1.3073 2.090E-008 
6.900E-001 1.3071 2.400E-008 
7.000E-001 1.3069 2.900E-008 
7.100E-001 1.3067 3.440E-008 
7.200E-001 1.3065 4.030E-008 
7.300E-001 1.3062 4.300E-008 
7.400፻-001 1.3060 4.920E-008 
7.500Е-001 1.3059 5.870E-008 
7.600E-001 1.3057 7.080E-008 
7.700E-001 1.3055 8.580E-008 
7.800E-001 1.3053 1.020E-007 
7.900E-001 1.3051 1.180E-007 
8.000E-001 1.3049 1.340E-007 
8.100E-001 1.3047 1.400E-007 
8.200E-001 1.3046 1.430E-007 
8.300E-001 1.3044 1.450E-007 
8.400E-001 1.3042 1.510E-007 
8.500E-001 1.3040 1.830E-007 
8.600E-001 1.3039 2.150E-007 
8.700E-001 1.3037 2.650E-007 
8.800E-001 1.3035 3.350E-007 
8.900E-001 1.3033 3.920E-007 
9.000E-001 1.3032 4.200E-007 
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X, мп n X 

9.100E-001 1.3030 4.440E-007 
9.200E-001 1.3028 4.740E-007 
9.300E-001 1.3027 5.110E-007 
9.400E-001 1.3025 5.530E-007 
9.500E-001 1.3023 6.020E-007 
9.600E-001 1.3022 7.550E-007 
9.700E-001 1.3020 9.260E-007 
9.800E-001 1.3019 1.120E-006 
9.900E-001 1.3017 1.330E-006 
1.000E + 000 1.3015 1.620E-006 
1.010E + 000 1.3014 2.000E-006 
1.020E + 000 1.3012 2.250E-006 
1.030E + 000 1.3010 2.330E-006 
1.040E + 000 1.3009 2.330E-006 
1.050E + 000 1.3007 2.170E-006 
1.060E + 000 1.3005 1.960E-006 
1.070E + 000 1.3003 1.810E-006 
1.080E + 000 1.3002 1.740E-006 
1.090E + 000 1.3000 1.730E-006 
1.100E + 000 1.2998 1.700E-006 
1.110E + 000 1.2997 1.760E-006 
1.120E + 000 1.2995 1.820E-006 
1.130E + 000 1.2993 2.040E-006 
1.140E + 000 1.2991 2.250E-006 
1.150E + 000 1.2990 2.290E-006 
1.160E + 000 1.2988 3.040E-006 
1.170E + 000 1.2986 3.840E-006 
1.180E + 000 1.2984 4.770E-006 
1.190E + 000 1.2982 5.760E-006 
1.200E + 000 1.2980 6.710E-006 
1.210E + 000 1.2979 8.660E-006 
1.220E + 000 1.2977 1.020E-005 
1.230E + 000 1.2975 1.130E-005 
1.240E + 000 1.2973 1.220E-005 
1.250E + 000 1.2971 1.290E-005 
1.260E + 000 1.2969 1.320E-005 
1.270E + 000 1.2967 1.350E-005 
1.280E + 000 1.2965 1.330E-005 
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Table A.1 (continued) 








































































































Appendix 

X, итп n X 

1.290E 4- 000 1.2963 1.320E-005 
1.300E 4- 000 1.2961 1.320E-005 
1.310E 4- 000 1.2959 1.310E-005 
1.320E 4- 000 1.2957 1.320E-005 
1.330E 4- 000 1.2955 1.320E-005 
1.340E 4- 000 1.2953 1.340E-005 
1.350E 4- 000 1.2951 1.390E-005 
1.360E 4- 000 1.2949 1.420E-005 
1.370E 4- 000 1.2946 1.480E-005 
1.380E 4- 000 1.2944 1.580E-005 
1.390E + 000 1.2941 1.740E-00 
1.400E + 000 1.2939 1.980E-005 
1.410E 4- 000 1.2937 3.442E-005 
1.420E + 000 1.2934 5.959E-005 
1.430E 4- 000 1.2931 1.028E-004 
1.440E + 000 1.2929 1.516E-004 
1.449E + 000 1.2927 2.030E-004 
1.460E 4- 000 1.2924 2.942E-004 
1.471E 4- 000 1.2921 3.987E-004 
1.481E + 000 1.2920 4.941E-004 
1.493E + 000 1.2918 5.532E-004 
1.504E + 000 1.2916 5.373E-004 
1.515E + 000 1.2914 5.143E-004 
1.527E + 000 1.2912 4.908E-004 
1.538E + 000 1.2909 4.594E-004 
1.563E + 000 1.2903 3.858E-004 
1.587E + 000 1.2897 3.105E-004 
1.613E 4- 000 1.2890 2.659E-004 
1.650E 4- 000 1.2879 2.361E-004 
1.680E 4- 000 1.2870 2.046E-004 
1.700E 4- 000 1.2863 1.875E-004 
1.730E 4- 000 1.2853 1.650E-004 
1.760E 4- 000 1.2843 1.522E-004 
1.800E 4- 000 1.2828 1.411E-004 
1.830E + 000 1.2816 1.302E-004 
1.840E + 000 1.2811 1.310E-004 
1.850E + 000 1.2807 1.339E-004 
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A, ym n X 

1.855E + 000 1.2805 1.377E-004 
1.860E + 000 1.2802 1.432E-004 
1.870E + 000 1.2797 1.632E-004 
1.890E 4- 000 1.2788 2.566E-004 
1.905E 4- 000 1.2780 4.081E-004 
1.923E 4- 000 1.2771 7.060E-004 
1.942E + 000 1.2762 1.108E-003 
1.961E + 000 1.2756 1.442E-003 
1.980E + 000 1.2750 1.614E-003 
2.000E + 000 1.2744 1.640E-003 
2.020E + 000 1.2736 1.566E-003 
2.041E + 000 1.2728 1.458E-003 
2.062E 4- 000 1.2718 1.267E-003 
2.083E 4- 000 1.2707 1.023E-003 
2.105E 4- 000 1.2694 7.586E-004 
2.130E 4- 000 1.2677 5.255E-004 
2.150E 4- 000 1.2663 4.025E-004 
2.170E 4- 000 1.2648 3.235E-004 
2.190E 4- 000 1.2633 2.707E-004 
2.220E + 000 1.2609 2.228E-004 
2.240E + 000 1.2591 2.037E-004 
2.245E + 000 1.2587 2.026E-004 
2.250E 4- 000 1.2582 2.035E-004 
2.260E 4- 000 1.2573 2.078E-004 
2.270E + 000 1.2564 2.171E-004 
2.290E + 000 1.2545 2.538E-004 
2.310E + 000 1.2525 3.138E-004 
2.330E + 000 1.2504 3.858E-004 
2.350E + 000 1.2482 4.591E-004 
2.370E + 000 1.2459 5.187E-004 
2.390E + 000 1.2435 5.605E-004 
2.410E 4- 000 1.2409 5.956E-004 
2.430E 4- 000 1.2382 6.259E-004 
2.460E 4- 000 1.2337 6.820E-004 
2.500E 4- 000 1.2270 7.530E-004 
2.520E + 000 1.2232 7.685E-004 
2.550E + 000 1.2169 7.647E-004 
2.565E + 000 1.2135 7.473Е-004 
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Appendix 

X, мп n X 

2.580E + 000 1.2097 7.392E-004 
2.590E + 000 1.2071 7.437Е-004 
2.600E + 000 1.2043 7.543Е-004 
2.620Е + 000 1.1983 8.059Е-004 
2.675Е + 000 1.1776 1.367Е-003 
2.725Е + 000 1.1507 3.508Е-003 
2.778E + 000 1.1083 1.346E-002 
2.817E + 000 1.0657 3.245E-002 
2.833E + 000 1.0453 4.572E-002 
2.849E + 000 1.0236 6.287E-002 
2.865E + 000 1.0001 8.548E-002 
2.882E + 000 0.9747 1.198E-001 
2.899E + 000 0.9563 1.690E-001 
2.915E + 000 0.9538 2.210E-001 
2.933E + 000 0.9678 2.760E-001 
2.950E + 000 0.9873 3.120E-001 
2.967E + 000 1.0026 3.470E-001 
2.985E + 000 1.0180 3.880E-001 
3.003E + 000 1.0390 4.380E-001 
3.021E + 000 1.0722 4.930E-001 
3.040E + 000 1.1259 5.540E-001 
3.058E + 000 1.2089 6.120E-001 
3.077E + 000 1.3215 6.250E-001 
3.096E + 000 1.4225 5.930E-001 
3.115E + 000 1.4933 5.390E-001 
3.135E 4- 000 1.5478 4.910E-001 
3.155E + 000 1.5970 4.380E-001 
3.175E 4- 000 1.6336 3.720E-001 
3.195E 4- 000 1.6477 3.000E-001 
3.215E + 000 1.6405 2.380E-001 
3.236E + 000 1.6248 1.930E-001 
3.257E + 000 1.6108 1.580E-001 
3.279E + 000 1.5905 1.210E-001 
3.300E + 000 1.5714 1.030E-001 
3.322E + 000 1.5559 8.360E-002 
3.345E + 000 1.5396 6.680E-002 
3.367E + 000 1.5241 5.312Е-002 
3.390E + 000 1.5086 4.286E-002 
3.413E + 000 1.4949 3.523E-002 
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Table A.1 (continued) እ, pm ü x 
3.436E 4- 000 1.4827 2.887E-002 
3.460E 4- 000 1.4710 2.347E-002 
3.484E + 000 1.4604 1.921E-002 
3.509E + 000 1.4502 1.586E-002 
3.534E + 000 1.4411 1.326E-002 
3.559E + 000 1.4328 1.130E-002 
3.624E + 000 1.4146 8.146E-003 
3.732E 4- 000 1.3924 6.672E-003 
3.775E 4- 000 1.3850 6.966E-003 
3.847E 4- 000 1.3750 8.248E-003 
3.969E 4- 000 1.3623 1.112E-002 
4.099E + 000 1.3526 1.471E-002 
4.239E + 000 1.3447 1.867E-002 
4.348E + 000 1.3406 2.411E-002 
4.387E + 000 1.3401 2.656E-002 
4.444E + 000 1.3412 2.990E-002 
4.505E + 000 1.3444 3.179E-002 
4.547E + 000 1.3473 3.090E-002 
4.560E + 000 1.3482 3.007E-002 
4.580E + 000 1.3491 2.883E-002 
4.719E + 000 1.3470 1.940E-002 
4.904E + 000 1.3379 1.347E-002 
5.000E + 000 1.3325 1.240E-002 





netration in snow and ice, estimation of actinic flux, etc. The results for the ice 
refractive index by different authors are given in Tables A.1 and A.2. The uncertainty 
of measurements can be understood from Table A.3 and Fig. A.1. 
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Table A.2 The spectral 
dependence of imaginary part 
of ice refractive index (Picard 
et al. 2016) 








































































































Appendix 

À, ит X 

0.320 8.3888e-10 
0.325 7.9573e-10 
0.330 7.9094e-10 
0.335 7.7136e-10 
0.340 7.47346-10 
0.345 7.14326-10 
0.350 6.99106-10 
0.355 6.78076-10 
0.360 6.63786-10 
0.365 6.55306-10 
0.370 6.47806-10 
0.375 6.39946-10 
0.380 6.42046-10 
0.385 6.27806-10 
0.390 6.33966-10 
0.395 6.16156-10 
0.400 6.27516-10 
0.405 6.06106-10 
0.410 5.84876-10 
0.415 5.71556-10 
0.420 5.76456-10 
0.425 5.86496-10 
0.430 5.97246-10 
0.435 6.19606-10 
0.440 6.37436-10 
0.445 6.61166-10 
0.450 6.89566-10 
0.455 7.18636-10 
0.460 7.5197e-10 
0.465 7.8996e-10 
0.470 8.4021e-10 
0.475 8.9213e-10 
0.480 9.4687e-10 
0.485 1.0089e-09 
0.490 1.0791e-09 
0.495 1.1593e-09 
0.500 1.2457e-09 
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Table A.2 (continued) እ, pm " 
0.505 1.3446e-09 
0.510 1.4568e-09 
0.515 1.5652e-09 
0.520 1.6960e-09 
0.525 1.8334e-09 
0.530 1.9857e-09 
0.535 2.1573e-09 
0.540 2.3470e-09 
0.545 2.5657e-09 
0.550 2.8008e-09 
0.555 3.0635e-09 
0.560 3.3461e-09 
0.565 3.6448e-09 
0.570 3.9598e-09 
0.575 4.3020e-09 
0.580 4.6519e-09 
0.585 5.0425e-09 
0.590 5.4761e-09 
0.595 5.9682e-09 
0.600 5.7300е-09 








Table EE ВА inue ኢ, wm | x (Warren 1984) | x (Warren апа | x (Picard et al. 
part of ice refractive index Brandt 2008) 2016) 


(multiplied by 10°) at 

















selected wavelengths in the 0.35 3.675 0.020 0.700 
visible from different sources. 0.40 2.710 0.024 0.628 
The differences in the spectral በ 45 1.540 0.092 0.690 
range 0.35-0.55 jum are 

particularly large x (Picard 0.50 1.910 0.589 1.246 
et al. 2016) 0.55 3.110 2.289 2.801 





0.60 5.730 5.730 5.730 
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Fig. A.1 Real (n) and imaginary (k) parts of ice (blue lines) and water (red lines) refractive indices 
(Gallet et al. 2014). There are several wavelengths, where the values of the imaginary parts of ice 
and water refractive indices coincide. One of them is located close at 1300 nm. The real part of ice 
refractive index is systematically smaller as compared to the value of п for water 


A.2 The Polarization Characteristics of Singly Scattered 
Light 


The Stokes vector of the reflected light beam 5 in the single scattering approximation 
can be presented in the following form (Hovenier et al. 2004): 


5 = M(x — фо) Р(Ө)М(—ф\)$у, (A.2.1) 


where Ë (0) is the 4 x 4 phase matrix, M (ф) is the 4 x 4 rotational matrix, So is the 
Stokes vector of the incident solar light. The rotational matrices in Eq. (A.2.1) are 
needed because the Stokes vectors So(S ) are defined with respect to the meridional 
plane containing the direction of light incidence (observation) and normal to the 
scattering layer. The phase matrix is defined with respect to the scattering plane 
containing directions of incident and scattered light. It follows for the rotational 
matrix: 


1 0 0 0 
0 cos 2¢@ — sin2¢ 0 
0 510 2ф cos2¢ 0 
0 0 0 1 


М(ф) = (A.2.2) 
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The angles фу and 2 can be derived from the following equations (Hovenier et al. 
2004): 


—Ho ተ Huy —H + Hoy 


, COS Q» = | 
sJ/ü- 430 — у?) sy — 490 - у?) 


where шо is the cosine of the incidence zenith angle, и is the cosine of the viewing 
zenith angle, s = sgn(w — л), v is the relative azimuthal angle and 


y = uuo + 0. — i — u?) sin y (ላ.2.4) 


is the cosine of scattering angle. The phase matrix has been introduced in Chapter 3. 
It follows at #0 = 1: 


cos фу = 








(A.2.3) 


y = u; COSY] = COS #2 = 1. (A.2.5) 


Then rotational matrix coincides with the unity matrix. Assuming that the incident 
light is unpolarized, one derives from Eq. (A.2.1): 


I Py lo 

Q Е cos(2¢) Pra fo | (А.2.6) 
U — sin(292) Pi Io 

V 0 


Therefore, it follows for the degree of linear polarization: 


p=,/p2 + pi, (A.2.7) 





where 
Q Р 
= = — cos 2 A.2.8 
Pq U p2 Р ( ) 
апа 
U . Pj 
u = — = — sin 2g ——. A.2.9 
p I #2 Р ( ) 


One can see that p, = tg(295) р „р = p, for the case considered. It follows that 
in the principal plane: #2 = 0 or л, tg (292) = 0 and, therefore, p; = — pi» = — Fr 


Pu 
and p, = 0. 
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А.З The Simplified Radiative Transfer Model 


We shall present the top-of-atmosphere (TOA) reflectance for atmosphere-underlying 
snow system in the following way (Liou 1992, 2002; Kokhanovsky 2020): 


Ta fr; 


1 — rar; 





Ктол = RT, = Ra + : (A.3.1) 


where R is atmospheric reflectance of an idealized gas—free terrestrial atmosphere, 
R, the atmospheric (with account for molecular scattering and aerosol absorption and 
scattering processes) contribution to the TOA reflectance, r, is the spherical albedo 
of the atmosphere, r, is the bottom-of-atmosphere snow spherical albedo, T, is the 
total atmospheric transmittance from the top-of-atmosphere to the underlying surface 
and back to the satellite position (without account for the gaseous components of the 
atmosphere), T, is the atmospheric transmittance due to the gaseous components of 
the atmosphere. 

In the case of Lambertian underlying surfaces, the underlying surface reflectance 
does not depend on solar and viewing observation directions and Eq. (A.3.1) is valid 
with f = / and r, = R;, where R, is underlying Lambertian surface reflectance. The 
snow is not exactly Lambertian reflector, therefore, we introduce the factor f in Eq. 
(A.3.1) to partially account for the non-Lambertian character of snow reflectance. 
We shall assume that 


Rsnow 
f2——, (A.3.2) 


Fs 


where Ко, is the snow reflection function. 

The top of atmosphere reflectance R, for a clean atmosphere can be presented in 
the following way using the Sobolev approximation derived from the scalar radiative 
transfer equation (RT) for non-absorbing media (Sobolev 1972; Kokhanovsky et al. 
2020): 


Ra = Rss + Rms, (A.3.3) 


where single scattering contribution 


Rss = М(т)р(0) (A.3.4) 


and multiple light scattering contribution is approximated as 


N(t) 


Rns =1+M i) M тавче 
+ M(t)q (Ho, и) 4430 — ge 


(A.3.5) 
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where 
1— e m 
— — = A. ` 
M(t) dup e fuo) (ш), (A.3.6) 
3 3 -t -1 -1 
f()-ltjut(l-gnJe*m-—nu, tH › (A.3.7) 
q (Mo, ш) = 3(1 + 8) MoM — 2000 + и). (А.3.8) 


Here, 4o is the cosine of the solar zenith angle (SZA), jz is the cosine of the 
viewing zenith angle (УХА), 0 is the scattering angle defined as cos0 = -ይዐይ ተ 
Sos COS 9, where g is the relative azimuthal angle (equalto 180 degrees minus OLCI 
relative azimuthal angle), so is the sine of the SZA, s is the sine of the УЛА, ፲ is 
the atmospheric optical thickness, p(0) is the phase function, g is the asymmetry 
parameter determined by the following expression: 


g= ; f p(0) sin Ө соѕ 040. (A.3.9) 
0 
The approximate account for aerosol absorption effects can be performed multi- 
plying Rss by the single scattering albedo wọ. The accuracy of Eqs. (А.3.1-А.3.3) 
can be further improved using the truncation approximation as discussed by Katsev 
et al. (2010). 
The phase function is presented in the following form (assuming that absorption 
effects can be ignored): 








p(0) == (Taer Paer (0) F Tn Pm(0))/T, (A.3.10) 
where 

Paer (0) = CPaer,1 (0) F (1 = C) Paer.2(), (A.3.11) 

1— TS 1= Baer 
Paeri (0) = - 3^ Paer2 (0) = 2 3^ 

(1 — 28aer,1 COS Ө + 821) ? 1 — 28aer,2 COSA + 81.2) ? 

(A.3.12) 
8аег,1 = 0.8, Saer2 = —0.45, (A.3.13) 
ር= M Baer? (А.3.14) 


, 
S8aer,1 — Saer,2 


px (0) = (1 + cos? 0) (A.3.15) 
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is the molecular scattering phase function. Therefore, it follows for the asymmetry 
parameter: 


Taer 
CEN a, A.3.16 
5 Tmol ተ Taer 5 ( 


The parameter g, varies with the location, time, aerosol, type, etc. We shall assume 
that it can be approximated by the following equation: 


Ba = 80 + me s. (A.3.17) 


The coefficients in this equation (as derived fro ultiple year AERONET 
observations over Greenland (Kokhanovsky et al. 202 е: 


go = 0.5263, g, = 0.4627, ào = 0.4685 177. (A.3.18) 
The transmission function T (uo, ш) is approximated as follows: 


T (uo, ш) = t", (A.3.19) 


where t is calculated using the following approximation (Katsev et al. 2010): 


t=e BT, (A.3.20) 


where 


f ዖ(ዐ) sinodo (A321) 


B= 


NI = 


л 
2 


is the so -- called backscattering fraction and т is the atmospheric optical thickness. 
It follows from Eqs. = and-A&3-10): 


В = (tae; Ba + Tm Bm)/T, (A.3.22) 


where B,, — 0.5, 


B, = C Baer (Sani) + (1 = C) Baer (gaer,2), (A.3.23) 


Baer (g) = x lrg ) (A.3.24) 


у1+ 82 
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Table A.3.1 The polynomial expansion coefficients of the parameters given in Eq. (A.3.26). The 


3 
parameter M is given by the following equation: М = ኔ Mgg* (and similar for other parameters) 

















s=0 
5 0 1 2 3 
M; 0.18016 —0.18229 0.15535 —0.14223 
Ns 0.58331 —0.50662 —0.09012 0.020700 
Ds 0.21475 —0.10000 0.13639 —0.21948 
Es 0.16775 —0.06969 0.08093 —0.08903 
Ks 1.09188 0.08994 0.49647 —0.75218 

















It follows from Eq. (A.3.24) that Baer (1) = 0 as it should be. Also we have from 
Eq. (A.3.24) at small values of g: 





1 2—3 
Baer (g) = + sls ) 


2 20+ 9? +(1— 8) V1+ 82). 


Therefore, it follows: Baer(0) = 0.5 as it should be for the symmetric phase 
functions (e.g., for molecular scattering). 

The atmospheric spherical albedo r, is found using the approximation proposed 
by Kokhanovsky et al. (2005): 


(A.3.25) 


yes (uet +Ne + D) т: (A.3.26) 


The coefficients of polynomial expansions of all coefficients (M, N, D, £, к) in 
Eq. (A.3.26) with respect to the value of the asymmetry parameter g are given in 
next Table A.3.1 (Kokhanovsky et al. 2005). 

One can see that the functions discussed above depend on the atmospheric optical 
thickness, which can be presented in the following form: 


T(A) = Tmol (A) + Taer (A). (A.3.27) 
The molecular optical thickness can be approximated as (Iqbal 1983): 


Umol (A) = qu” (A.3.28) 


at the normal pressure po and temperature To. Here, q = 0.008735, v = 4.08, and 
the wavelength is in microns. The value of molecular optical thickness at another 
pressure level p can be derived using the following expression: Tmo (À) = PT (A), 
where p — 5 , p is the site pressure, ро = 1013.25mb. The site pressure is calculated 


as: p — po exp( — £) . Here z is the height of the underlying surface and H — 7.64 km 
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is the scale height. The molecular optical thickness at Arctic and Antarctic sites is 
discussed in detail by Tomasi and Petkov (2015). 
The aerosol optical thickness can be presented as 


À 


Taer (À) = (=) , (A.3.29) 


where Ao = 0.4um, (œ, В = Taer (Ao)) are the Angstróm parameters. One can assume 
that о = 1.699% et al. 2005). 

As far as us transmission is of concern, one can use the following exponential 
approximation: 


Т,(А) = exp(—At,(A)), (А.3.30) 


where А is the air mass factor (Rozanov and Rozanov 2010) апа т, (A) is the gaseous 
absorption optical thickness. 

The air mass factor can be approximated by the following equation for the case 
of ozone (Iqbal 1983): 


А 1+5 1+5 
J2s ud УО +p? 





(А.3.31) 





s — H/a, Н is the height of ozone layer (assumed to be equal 22 km), а = 6370 km 
is the radius of the Earth. The account for the latitudinal change of the ozone layer 
height can be performed using the following empirical relationship (Savastiouk and 
McErloy 2004): H = x — 0.1y, where y is the latitude in degrees without regard for 
the sign (e.g., 90? for the South and North Poles) and x — 26 km. 
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